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TECHNICAL MEMORANDUM X-64781 


A COMPARISON OF DIGITAL COMPUTER PROGRAMS FOR THE 
NUMERICAL SOLUTION OF ORDINARY 
DIFFERENTIAL EQUATIONS 


When a physical problem is to be simulated with the aid of a digital 
computer, the result is often a system of ordinary differential equations that 
must be solved numerically. For this reason an efficient numerical integra- 
tion algorithm is very desirable to minimize the amount of computer time 
needed to solve a particular problem. 

Many different approaches are available for the numerical solution of 
system of ordinary differential equations, and an evaluation in this report of 
some of these approaches was motivated by a desire to compute optimal 
trajectories as rapidly as possible. In computing optimal trajectories, the 
system of ordinary differential equations known as the equations of motion and 
adjoint variables must be integrated many times to satisfy iteratively the 
boundary conditions. This process consumes much computer time unless the 
particular numerical integration technique being used is operating efficiently. 

The Runge-Kutta formulas derived by E. H. Fehlberg in References 1 and 2 
and discussed in the next sections of this report meet this requirement of 
operational efficiency better than any of the other methods tested. 

Also described in the section of this report entitled An Explanation of 
Runge-Kutta Numerical Integration Formulas is a new stepsize control proce- 
dure for Fehlberg's numerical integration formula. This new stepsize control 
procedure is combined with one of the best of Fehlberg’s formulas, and the 
results it achieves for a variety of test problems is compared with the results 
achieved by other numerical integration techniques for the same test problems. 
The data obtained from these comparisons is discussed in detail in the sections 
of this report entitled Fehlberg’s Runge-Kutta Formulas With Stepsize Control 
and A Comparison of Fehlberg’s 7-8-13 Formulas With the Numerical Integration 
Techniques of References 3 and 4. 



AN EXPLANATION OF RUNGE-KUTTA NUMERICAL 
INTEGRATION FORMULAS 


Consider the system of differential equations denoted by 


x = f (t, x) , (l) 

where x, f, and x are all vectors of dimension n. Note that the symbol x 
denotes dx/dt where t is assumed to be the independent variable for the sys- 
tem. For a given set of initial conditions denoted by x (t Q ) = values for 

x(t Q + At) can be obtained from Runge-Kutta formulas as follows: 


m 


x(t o+ At) = x o + At £ , 

k=0 


( 2 ) 


where 


f = f (t , x ) 


o o 


k-1 


£r = f (t 0 + a k At, x o + At I t t ) for k= 1, 2, . . . ,m . (3) 


The coefficients a^, c^, /3^ , and the integer m are determined to make the 

expression for x(t Q + At) as given in equation (2) equal to a Taylor series 

expansion for x(t Q + At) up to a certain order. For this determination of the 

coefficients to be meaningful, the vector function f(t, x) must be reasonable 

enough to have a convergent Taylor series expansion in some neighborhood of 

the initial conditions t , x . 

o o 

As an example of how the coefficients c^, and (3^ are determined, 

a second-order Runge-Kutta formula can be derived. To do this, consider a 
Taylor series expansion of x(t Q + At) to second order. That is. 


2 



( 4 ) 


x (t + At) = x + x At + Vo x At 2 
o o o 2 o 

Then note that 

x = f (t , x ) 
o o o 

and 



o o 


Thus the Taylor series expansion can be rewritten as: 


x(t Q +At) X 0 + f (V X o) At+V * (dt),, + (fx) fJ ' 11 




At 


(5) 


Now the expression for x(t + At) from the Runge-Kutta approach is given by 
equation (2). That is, ° 


x (t + At) = x +At[c f (t , x ) + c, f. + ... + c f ] 
o o o o o 11 mm 


( 6 ) 


Now equation (3) can be used and fj expanded in a multivariable Taylor 
series to give: 


I: 


3 


*1 “ f (t o + “l ^ x o + ^loV 

- f <V x o> + (l) 0 (“l * (l) 0 [4t % V + - < 7 > 


When the above expression is substituted into equation (6) , second-order terms 
will be obtained in At. Thus, there is no need to carry the expansion of 
equation (7) past first order and m in equation (6) is chosen to be one. Thus 
equation (6) becomes: 


x (t Q + At) 


x + At { e f(t,x) 
o 1 o ' o o' 


+ c. ["f (t , x ) + (a. At) + (At/3, ft f )1 

1[^ v o o \8t/ 0 v 1 ' \d x/ 0 10 O' I 


( 8 ) 


Therefore, 


x (t + At) = x + (c + c,) f (t , x ) At 
v o ’ o v 0 1 v o o' 

+ [(If) K a i) + (tt) f (t , x ) c. fl 
|\ 8 t / 0 1 l ' \ 9x / 0 0 0 1 10 


At" (9) 


Now a comparison of equations (5) and (9) shows that: 


C 0 + C 1 = 1 


C i a i - '* 


= V, 


Vio = V * 


do) 


4 
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Any set of coefficients that satisfies equation (10) will give a second-order 
Runge-Kutta formula that uses only two evaluations of the differential equations. 
Thus it can be seen that Runge-Kutta formulas are not unique. As an example 
c 0 = Cj = V 2 , = 1, andj3 10 = 1 will satisfy equation (lo) . Equation (10) is 

called the equations of condition for a second-order Runge-Kutta formula. 

Also, the number of evaluations of the differential equations required for a 
particular order is usually called the number of function evaluations needed. 
With this background information, the Runge-Kutta formulas developed by 
Fehlberg in References 1 and 2 can now be discussed. 


FEHLBERG'S RUNGE-KUTTA FORMULAS 
WITH STEPSIZE CONTROL 


To obtain a useful stepsize control procedure for Runge-Kutta formulas 
some indication of the truncation error of the series expansion must be deter- 
mined. Fehlberg's idea is to develop Runge-Kutta formulas of adjacent order 
that use the same function evaluations, and then the difference in the two 
formulas is a good approximation to some single term in the Taylor series 
expansion. 


In References 1 and 2, Fehlberg performs the very difficult determina- 
tion of adjacent Runge-Kutta formulas for orders from one to eight. That is, 
formulas of order 1 and 2, 2 and 3, 3 and 4, 4 and 5, 5 and 6, 6 and 7, 7 and 
8, and 8 and 9 are all developed, and each adjacent pair uses the same function 
evaluations so that their difference can approximate a corresponding term in 
the Taylor series expansion for x(t + At). As an example, for m= 12 a 
seventh-order Runge-Kutta formula and an eighth-order Runge-Kutta formula 
are developed that use the same function evaluations. That is. 


12 

x(t +At) = x + At 2 c f (11) 

o o i^o 


x (t + At) 
o 


X + At 
o 



c k 


( 12 ) 


5 



where 


f 0 = f <«o> x o> 


k-1 


f k = f (t Q + o ' k At, x q + At Yj ^ f *) for k = 1, 2 12 . (13) 


Table 1, which is taken from Reference 1, gives the values for c , 0 , at , and 

K K K 

/3 . Equation (ll) is a Runge-Kutta formula that agrees with the Taylor 

k* 

series expansion to order seven, and equation (12) is a Runge-Kutta formula 
that agrees with the Taylor series expansion to order eight. Thus, the differ- 
ence between the two formulas that Fehlberg denotes by TE is a good approxi- 
mation to the eighth-order term in the Taylor series expansion. This difference 
is given by the following expression. 


TE = 


41 

840 


(f + f - f - f ) At 

v 0 10 11 12' 


(14) 


Since the expression for TE as given by equation ( 14) is a good approximation 
to the eighth -order term in the Taylor series expansion for x(t Q + At) , it can 

be used to determine At. Let I TE| denote the largest component of the vector 
TE as given in equation (14). Then 


I TE I « a At 8 

O 


(15) 


where a 8 is the coefficient of the eighth-order term in the Taylor series ex- 
pansion for the variable which yields the largest component of TE. Therefore, 


a « 
8 


1 TE | 

At 8 


(16) 
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Thus, for a given value of At, equation (14) can be used to compute TE. Then 
equation (16) can be used to compute a 8 . Let e denote the largest acceptable 
component of TE . Then, if I TE I ^ e , the step taken with At is an acceptable 
step and a At for the next step can be computed. If I TE| > e, then the step 
taken with At is not acceptable, and the step must be taken over again with a 
new At thatf can be computed. Let At^ denote the new At that is to be com- 
puted and At Q denote the old At. To derive an expression for At^, it is 
required that 


e 


a At 
8 n 


(17) 


That is, the new At should produce a maximum component of the truncation 
error of magnitude e. Thus 



Now the value of a 8 is obtained from equation (16) with the value of I TE| 
produced by At^. Thus, 



Equation (19) gives the computed value of At^ that is used for the next integra- 
tion step if | TEl ^ e and is used for a repeat step if I TE| > e. In actual 
practice, only 80 percent of the At Q as computed by equation (19) is used for 

either the next step or a repeat step. Also, in computing | TEl each component 
is usually divided by the corresponding component from the initial value of x so 


8 



I 


that a valid comparison can be made in search of the maximum component. 
Several example computer listings are shown in Appendix A so that the different 
implementations of the normalization of the truncation error vector can be 
illustrated. 

Note that in the example computer program listings, equation (12) is 
used to compute x(t Q + At) . This allows more accuracy in the final result 

to be obtained than if equation ( ll) were used because equation ( 12) is an 
eighth -order formula. In fact, for most problems, e will be a bound for the 
entire solution time interval because of the conservative stepsize control pro- 
cedure used. Other less conservative procedures could be implemented if 
computer speed is of more importance than accuracy of the solutions obtained, 
but the results in the next sections indicate that this conservative stepsize 
control gives very satisfactory computer execution times on the example 
problems. The formula chosen for testing on the example problem is the 7-8 
formula which used 13 function evaluations (denoted by Fehlberg’s 7-8-13 
formula). Fehlberg states in Reference 2 that this formula is probably the 
best for difficult problems, but for easier problems where values of x(t) are 
needed at more frequent intervals, a lower order formula might be more 
efficient. If a lower order formula is used (for example, a 2-3 formula) , then 
equation (19) can still be used to compute At, but the power (l/8) is replaced 
by (1/3). 


A COMPARISON OF FEHLBERG'S 7-8-13 FORMULA 
WITH THE NUMERICAL INTEGRATION TECHNIQUES 
OF REFERENCES 3 AND 4 


From Reference 3, two example problems are selected. These 
examples are denoted by B1 and FI. The definition of problem B1 is given 
as follows: 


Test Problem B1 


X 2 X 
1 X 2 

X 1 <‘o> = 1 

t = 0 
o 


x 2 (t o' * 1 

ii 


9 


I 


The definition of problem FI is given as follows: 


Test Problem FI 


x„ = x 
1 3 


x„ = x„ 
2 4 


= X x + 2 x 4 - ( 1 - n) 


(x x +M) 


% 


'( \ 2 a. 2 1 /2 

(x 1 + **) + x 2 


- (m) 


(x x - 1 + M) 


[( x 1- 1 + “) 2+x 2 2 ] ’ 


X 4 = x 2 - 2 * 3 " ( 1 ~ 


- (m) 


r n3/ xrv r n3/ 

o o '2 o o ' 2 

(x x + m) +x 2 J [(Xj-I + m) + x 2 


x (t ) =0.994 
1 o 

x 2 (* 0 ) = 0 . 

X 3 ( V = °- 

x ( t ) = -2.0317326296 
4 v o 7 


M = 


1 + 

0.012277471 


t = 0 


t„ = 11.124340337 


Figure 1 shows the results obtained for test problem B1 using Fehlberg’s 
7-8-13 formula (denoted by RK713) . Figure 2 shows the results obtained for 
test problem FI using RK713. In Figure 1, the error shown is calculated as 
e = | e 4 - x (l) | . Also, note that the percentage error in x(2) is equal to 
the percentage error in x(l) because of the stepsize control for RK713 that 
was used for this problem. This can be seen by examining the computer pro- 
gram listing shown in Appendix A and the associated output for this problem. 

Figure 5 of Reference 3 is comparable to Figure 1 of this report. That 
is. Figure 5 of Reference 3 shows the same type of information for the integra- 
tion techniques tested in Reference 3, that is shown by Figure 1 of this report 
for RK713. From Figure 5 of Reference 3, it can be seen that a predictor- 
corrector technique (denoted by HPCG) uses the fewest function evaluations 
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100 


1 l i I 

200 300 400 500 

FUNCTION EVALUATIONS 


Figure 1. Test problem Bl. 





(about lOO) for an e of 10' 1 . As seen in Figure 1, RK713 needs about 85 
function evaluations for this accuracy. For an accuracy of e = 10“ 2 , RK713 
needs only 120 function evaluations, whereas the only two integration techniques 
tested in Reference 3 that were able to achieve this accuracy (rational extrap- 
olation and extrapolated Runge-Kutta) needed about 300 function evaluations. 

As seen in Figure 1, RK713 is able to achieve an e of 10 -4 with about 240 
function evaluations. None of the integration techniques tested in Reference 3 
was able to achieve this accuracy. The fact that the curve levels off at around 
10" 5 accuracy is due to the precision limits of the computer used and not the 
difficulty of the problem or any inadequacy in the RK713 technique. The 
computer runs in this section were made on an XDS-930 computer unless other- 
wise indicated. This computer is able to carry only 11 or 12 significant 
figures for each computation. The dotted line in Figure 1 shows the results 
from the UNIVAC 1108 computer using double precision, which allows 16 to 18 
significant figures to be considered in the computations. 

When Figure 2 of this report is compared with Figure 6 of Reference 3, 
the superiority of RK713 is again demonstrated. The error e shown in Figure 
2 of this report and Figure 6 of Reference 3 is computed by the formula 


e 




(t o ) ] 2 +[x 2 ( t f )] 



RK713 is able to attain an e = 10 -6 with 1600 function evaluations. Most of the 
integration techniques tested in Reference 1 could not achieve an accuracy of 
even e = 10~ 5 . The only two routines to come close were again the extrapolated 
Runge-Kutta routine and the rational extrapolation routine, and both required 
about 2500 function evaluations to get close to the accuracy e = 10 -5 . Again a 
sample computer program listing and its output are included in Appendix A for 
this problem. 

From Reference 4, two example problems were also selected for com- 
parison purposes and denoted by the symbols B12 and E22. The definition of 
problem B12 is given as follows: 

Test Problem B12 (this is also problem E of Reference 3) 


2 (X ! - X 1 X 2 ) 

X 1 ' 1 

t = 0 
o 

- X 2 (1 -V 

* 3 

t f = 20 
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The definition of problem E22 is given as follows: 


Test Problem E-22 



<‘o> 


= 2 


x = (1 -X, 2 ) x o - x x (t ) = 0 

2 1 ’ 2 1 2 v o' 


t = 0 
o 

t. = 20 
f 


Figures 3 and 4 are plots of the errors achieved by RK713 versus the number 
of function evaluations needed for problems B12 and E22, respectively. To 
attempt to compare with the results shown in Reference 4, the error e shown 
in Figures 3 and 4 is computed as follows: 


max. [y 1 (t f ) - (t f ), y g (t f ) - x g (t f )] 


where [y (t ), y (t )] are assumed to be the actual solution of the problems 

XX Ci I 

and are obtained by a very accurate integration of the two example sets of 
differential equations on the UNIVAC 1108 using double precision arithmetic. 

From Figure 3 it can be seen that RK713 on problem B12 needs 600 
function evaluations to achieve an accuracy of e = 10~ 3 , 915 function evaluations 
to achieve an accuracy of e = 10 -6 , and 1615 function evaluations to achieve an 
accuracy of e = 10~ 9 . From Reference 4, the best method tested (with respect 
to the actual machine time needed to solve problem B12) is the Bulirsch-Stoer 
method. The Bulirsch-Stoer method needed 992 function evaluations for e = 10 -3 , 
1873 function evaluations for e = 10 -6 , and 3128 function evaluations for e = 10 -9 . 
Since the computer program DETEST mentioned in Reference 4 was not avail- 
able, the other statistics given in Reference 4 were not compared. It is hoped 
that the authors of Reference 4 will be interested enough in RK713 to try it on 
their DETEST program. 

From Figure 4, it can be seen that RK713 on problem E22 needs 665 
function evaluations to achieve an accuracy of e = 10 -3 , 915 function evaluations 
to achieve an accuracy of e = 10 -6 , and 1625 function evaluations to achieve 
an accuracy of e = 10 -9 . Again, the best method tested in Reference 4 (with 
respect to the actual machine time needed to solve problem E22) was the 
Bulirsch-Stoer method. On this problem, the Bulirsch-Stoer method needed 
1109 function evaluations for e = 10 -3 , 2009 function evaluations for e = 10 -6 , 
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and 3588 function evaluations for e = 10 -9 . Note again that a sample computer 
program listing and its output for both of the examples from Reference 4 are 
also contained in Appendix A. For each of the example problems of this 
section, the last output from each example computer program listing shown in 
Appendix A is a detailed print so that the actual behavior of each function can be 
seen. This detailed print was integrated as accurately as possible on the XDS 
930 computer and also serves to illustrate the detailed print feature of the 
computer program. 

All of the results of this section indicate that RK713 is superior to any 
of the other methods tested in References 3 and 4. Since all the errors in this 
section were absolute errors from the actual solutions of each of the differential 
equations, a more meaningful comparison is attempted in the next section. 

That is, RK713 and a rational extrapolation algorithm DIFSYF are compared 
on the basis of the accuracy each integration algorithm thinks it is achieving. 
Since actual accurate solutions of most differential equations attempted are 
not available, this is really the only comparison that can be made for these 
problems. Another reason for the inclusion of the comparison of the rational 
extrapolation algorithm with RK713 is that References 3 and 4 both concluded 
that the rational extrapolation algorithm as developed by Bulirsch-Stoer in 
Reference 5 was probably the best all-around algorithm they tested. 


MORE DETAILED COMPARISONS OF FEHLBERG'S 7-8-13 
FORMULA WITH AN EXTRAPOLATION ALGORITHM 
BASED ON RATIONAL FUNCTIONS 


In this section, a comparison of the two integration techniques (Fehlberg’s 
7-8-13 formula RK713 and the rational function extrapolation algorithm DIFSYF) 
is made using systems of differential equations describing the following example 
problems: planar elliptic vacuum orbits with eccentricities ranging from 0.001 
to 0.991, a three-dimensional almost circular vacuum orbit, and two reentry 
trajectories with aerodynamic forces and the adjoint differential equations 
included. 1 

After preparation of this report was complete the author was informed 
that a newer version of DIFSYF had been developed by R. Bulirsch. This new 


1. The comparisons shown were performed by Dr. E. D. Dickmanns while 
he was working at Marshall Space Flight Center. They are included in this 
report with his permission. Dr. Dickmanns is now working in West Germany 
at the DFVLR - Institut fur Dynamik der Flugsysteme. 
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version is claimed to be 30 to 40 percent faster than the old version for some 
accuracy requirements. Appendix B contains a listing of both the old and new 
versions of DIFSYF. No comparisons were made with the new version of 
DIFSYF because the results shown in this section and the previous section in- 
dicate that a 30- to 40-percent improvement in DIFSYF would not make it 
superior to RK713 for most accuracy requirements. 

The following sections describe in detail the example problems and the 
comparisons that were made between RK713 and the old version of DIFSYF. 


Systems of Differential Equations Used 

Planar Elliptic Orbits. For the range angle from pericenter 9 (eccen- 
tric anomaly) as independent variable the equations of motion are 



dr _ GM 
d0 " ur 


dr _ rr 
d0 ~ u" 


u = horizontal velocity 

r = radial velocity 

r = radial distance from center of 
gravity 

GM *= 3.986032* 10 14 (metric units) 

Earth mass times gravity constant. 


Three-Dimensional Vacuum Tra jectories. In this case time has been 
chosen as independent variable. The equations then are 


du _ ur 
dt " " r 


dx 

dt 

dr 

dt 

dr 

dF 


= — sin x tan X 
r 


u 2 GM 



= r 


X = azimuth from north (positive east) 
X = geocentric latitude 
A = longitude 

co = rotational speed of the earth 
e 


18 



dA 

dt 


u 

— cos x 
r 


dA _ u sin x 
dt r cos A ~ e 


Reentry Equations for State and Adj oint Variables. These equations 
were written in a flight path-oriented axis system for a spherical, nonrotating 
earth and an exponential density-altitude relationship. 


Sp o 

v ■= -g Star - (C D +kc L “) v 2 e 


-f>h 


X = - 


C T Sp fl , . 

v , . L O -ph sin p 

— cos y cos x tan A + — ve 

r 2m cos y 


e 


v cos y 
r cos A 


cos X 


A 


v 

— cos y sin X 


h = v sin y 


r = R E + h 


The differential equations for the adjoint variables are rather lengthy 
and are not reproduced here, but are shown in Reference 6. The controls 
C and p are determined from 

Xj 


sin p = ; os = [(A Y /cosy) 2 + A 2 ] 

03 cos y * y 

C = [ —o) / ( vA kn)] l/(n_l) 

L v 
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These are highly nonlinear equations sensitive to small changes in 
some of the adjoint variables. 

Symbols used: 



k 

m 


n 

R 


E 


S 


v 

J3 


7 

e 

i 

M 

P o 

X 


Zero-lift drag coefficient (set to 0.04) 

Lift coefficient (control) 

Local gravitational acceleration 

Altitude above sea level 

Factor for lift-dependent drag (0.8 and l.o) 

Vehicle mass 

Power for lift-dependent drag (1.75) 

Earth radius (6371 km) 

g 

Reference area (— = 100 [kg/m 2 ]) 
m 

Inertial velocity 

Inverse density scale height (assumed (3= 1.45- 10 
tm' 1 }) 

Flight path angle 

Range angle in initial flight direction 
Adjoint variable to i 
Bank angle (control) 

Density for h= 0 (assumed 1.54 [kg/m 2 ]) 

Azimuth relative to initial flight direction 
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Cases Compared 

To investigate the stepsize adjustment properties of the integration 
routines, system 1 has been integrated for eccentricities e from 0.001 (almost 
circular) to 0.991 (highly elliptic) with the following initial conditions: 


u (°) = t~“(l + e )l [m/s] 
r(o) = 0 

r(o) = 6 571000 [m] 


Two complete orbits have been integrated for each set of initial conditions and 
accuracy requirement. 

For the system of differential equations describing three-dimensional 
motion in vacuum about a rotating earth, an almost circular orbit inclined by 
30 deg to the equatorial plane has been chosen: 


u = 7850 [m/s] 
o 

X q = 90 deg 

r = 150 [m/s] 
o 

r = 6 470 000 [m] 
o 

X = 30 deg 

A = 0 
o 

Four complete revolutions have been integrated with the period P 
given by 


II-GM 

n/ 2* E 3 


E = 


(u 2 + r 2 ) 
o o' 
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With the reentry equations a single three-dimensional skip has been 
chosen as test case. The initial conditions were 


V = 

o 

7950 

Av = -0.543197 
o 

o 

II 

0 

Ax = -2054.06 
o 

11 

o 

3 deg 

Ay = -46.1287 
o 

A = 
o 

0 

AA = 538.79 
o 

1! 

o 

A 

95 000 

Ah = -0. 000703044 
o 


Discussion of Results 

For weakly elliptic planar orbits (e ~ 0.25) DIFSYF requires about 
30 to 40 percent more evaluations of the differential equations than RK713. 

The average stepsize is about five times as large for DIFSYF as for RK713 
(Figs. 5 and 6) . 

For accuracies (tolerance per step) required of 10 -7 , one complete 
revolution is integrated in two steps by DIFSYF. As the eccentricity increases 
the optimal stepsizes for integrating the apse-portions and the flanks of the 
orbit become more and more different. For these cases, the smaller stepsize 
of RK713 allows a faster stepsize adjustment, resulting in increasing superi- 
ority of RK713 (Figs. 7 and 8). For highly eccentric orbits, RK713 requires 
less than half the evaluations of the right-hand side of the differential equations. 
The ratio of maximum to minimum stepsize is about equal for both methods 
(~ 15 for e = 0. 991) . RK713 seems to be more tolerant against too low accu- 
racy requirements. To obtain approximately periodic orbits, tolerances re- 
quired for RK713 were ~ 10~ 3 while DIFSYF required 10 -5 (Fig. 9) . 

For the almost circular three-dimensional vacuum orbit DIFSYF takes 
two integration steps per revolution for 1CT 5 tolerance per step compared to 
10 -7 for planar orbits. Here the sinus-like time traces for x and A require 
stepsize adjustment, which causes RK713 to become more superior as the 
accuracy requirement is increased. The average stepsize for DIFSYF is four 
times that of RK713. For high accuracies DIFSYF requires twice the number 
of function evaluations as RK713 (Fig. 10) . 
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TOLERANCE/STEP TOLERANCE/STEP 


Figure 7. Performance comparison RK713-DIFSYF, planar orbits — e= 0.496 
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Figure 8. Performance comparison RK7 13 -DIFSYF, planar orbits — e= 0.7435 
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Figure 10. Three-dimensional vacuum orbit (almost circular) . 





The picture is very similar for the integration of the reentry equations. 
These test runs were made on a UNIVAC 1108 in double precision for accura- 
cies per step of 10 -14 to 10" 2 . For a three-dimensional skip (45-deg plane 
change) the results are shown in Figure 11. In the usually taken accuracy range 
10 -5 to 10 -1 °, RK713 is slightly superior to DIFSYF. For high accuracy re- 
quirements, this superiority becomes more pronounced. It again seems to be 
a question of stepsize adjustment. 

To test the properties of the integration routines with a very demanding 
trajectory, k has been changed from 1.0 to 0.8, everything else remaining 
constant. 

The resulting trajectory is a highly oscillating spiral dive with the 
flight path angle varying from -2 to -55 deg and normal accelerations of 150 g. 
The results are shown in Figure 12. Again RK713 is superior and takes about 
one-fourth the stepsize of DIFSYF. It also shows a more benevolent behavior 
against low accuracy requirements for this numerically unstable system of 
differential equations. DIFSYF breaks down for accuracies lower than 10 -9 
while RK713 yields good results up to 10 -4 . The breakdown 2 seems to be re- 
lated to the minimum stepsize used (lower part of Figure 12). Also, the 
conservative stepsize control for RK713 discussed in the section of this report 
entitled Fehlberg’s Runge-Kutta Formulas With Stepsize Control contributes 
to the more stable behavior of RK713 for low accuracy specifications. 

Check computations with simple precision (eight digits) for tolerances 
greater than 10 -7 were made with RK713. Supposedly because of roundoff 
errors the number of function evaluations was slightly higher, while the mini- 
mum stepsize was slightly smaller. Breakdown occurred at the same tolerance 
requirement of 10“ 3 . 


2. Breakdown means that when evaluating the right-hand side of the differential 

equations numbers larger than allowed for in the standard library functions 
occurred. 
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NUMBER OF STEPS TAKEN NUMBER OF RHS EVALUATIONS 

FOR 150-SEC FLIGHT TIME FOR 150-SEC FLIGHT TIME 



— ► TOLERANCE/STEP 


x 

/ 



— ► TOLERANCE/STEP 


Figure 11. Comparison RK7 13 -DIFSYF, reentry equations — 
three-dimensional skip trajectory. 



NUMBER OF RHS EVALUATIONS FOR 

MINIMUM STEPSIZE (SEC) 150-SEC FLIGHT TIME 



10" 2 10' 3 


SIMPLE 
DOUBLE PRECISION 


_L 


_L 


10 


n -e 


10 " 10 " 10 ' 10 10 “ 
TOLERANCE/STEP 


10‘ 10 10' 11 


10 


12 


10 


-13 



Figure 12. Comparison RK713-DIFSYF, reentry equations — 
oscillatory spiral dive [highly demanding trajectory (K= 0.8)1. 

31 




CONCLUSIONS 


The results of this report indicate that Fehlberg’s 7-8-13 (RK713) 
formula with the stepsize control developed in this report is superior to all of 
the other techniques with which it was compared. That is, it was able to solve 
all the example problems as fast or faster than any of the other techniques for 
the entire range of integration accuracies tested. For higher accuracies the 
superiority of RK713 is particularly evident and in fact RK713 was able to 
achieve better accuracies on all of the problems than any of the other techniques 
with which it was compared. Since the rapid solution of systems of differential 
equations is extremely important to the author of this report, he would 
appreciate being informed by any of the readers of more efficient techniques 
for solving these problems. 
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APPENDIX A 


SAMPLE COMPUTER PROGRAM LISTINGS AND SAMPLE 
OUTPUT USING THE FEHLBERG 7-8-13 FORMULA 

The following list of symbols is given in order to facilitate the use of 
the programs: 


DTGI 

An estimate of the initial t to be used by the program. 
After this DTGI is used once, the program is able to 
compute its own At so that the value of DTGI used has 
practically no effect on the number of function evalua- 
tions needed for a solution of a particular problem. 

DTPI 

The delta print times for the full print option. That is, 
if FPTI < TF, then intermediate prints will occur at 
intervals of DTPI after the print at FPTI. 

FPTI 

The control for the full-print option of the program. 

If FPTI ^ TF, then only the initial and final time are 
printed. If FPTI < TF, then FPTI is the time of the 
first intermediate print. 

KT 

An input control number that allows the program to 
print every KT integration steps if FPTI and DTPI are 
both larger than TF. If FPTI and DTPI are less than 
TF, then KT must be a large number. 

TF 

The final time. 

TI 

The initial time. 
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TABLE A-l. TEST PROBLEM B1 

AASSIGN 

S»f-T0i 

SI CR/B8 MT1/L0 LP. 

aPEWIND 

MT1 • 


AF0RTRAN 

30/ |_fi 

• 


1 


DIMENSION X[23 


?. 


DIMENSION ALPH[133/3ETAri3/123/CH[133 


3 


COMMON ALPH/ BETA# CH 


4 

20 

READ 11/ TIjDTGI#T0L 


5 


READ 11/XC13/XC23 


6 


READ 11/TF/FPT1/DT D 1 


7 


READ 2C25/<T 


8 

2025 

FORMAT CI43 


9 

C 

constants for integration subroutine 


10 


DO 60 1*1/13 


11 


DO 50 J= 1 / 12 


12 

50 

BETA [ I / J] *0. 


13 


ALPH CID *0. 


14 

60 

cum =o. 


15 


CH r 6] =34,/105. 


16 


CH T73 =9./35. 


17 


CH [83 *CH [71 


18 


CH [93 * 9 • /280 • 


19 


CH C 10] =CH C9] 


20 


CH 112] *41 ./840. 


21 


CH Cl 3] *CH [12] 


22 


ALPH C23 *2./27. 


23 


AL°H C31 *1 • /9 • 


24 


ALPH C4] » 1 • /6 • 


P5 


Al°H [5] *5./l2. 


26 


ALPH C6] * .5 


27 


AL d H [73 *5./6* 


28 


ALPH [8] * 1 »/6 • 


29 


AL D H C93 *2 */3 • 


30 


ALPH CIO] *1 //3. 

3 

31 


AL n H C 1 1 ] *1. 


32 


ALPHC13]«1. 


33 


BETA C2/ 1 ] 8 2 « /27 ♦ 


34 


BETA [3/ 13*1 1 /36 • 


35 


BETA [4/ 13 *l./24t 


36 


BETA [5/ 13 *5./12. 


37 


BETA C6/ 13 * .05 


38 


BETA [7/13 »-25./108. 


39 


BETA C8/ 1 3 «31 »/300» 


40 


3ETA C9/ 13 =2. 


4 1 


BETA CIO/ 13 «-91./108t 


42 


BETACU/ 13 =2383./4100. 


43 


BETA [12/ 13 =3./205. 


44 


BETA [13/ 13 8 - 1 777* /4100 • 


45 


BETA [3/ 23 8 1 • /1 2 • 


46 


beta [4/33 «l./8. 


47 


BETA C5/33 8 -25./16. 


48 


BETA [5/ 43 8 -BETA[5/33 


49 


beta [6/43 8 • 25 . 


rQ 


BETA C7/ 43 8 1 25 • /108 • 
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TABLE A-l. (Continued) 


51 


SETA [9/43 =-53. /6. 

52 


SETA [10/43 =23./108. 

^3 


BETA[ll/4]=-341./164. 

54 


BETA [13/4] =BETA Cll/43 

55 


BETA [6/53 « • 2 

56 


BETA [7/53 =-65. /27. 

57 


BETA [8/53 =61 ./225. 

58 


BE T A [9/53 =704. /45. 

^9 


BETA [10,53 »-976./135. 

6 0 


BETA [11/53 = 4 496./ 1025. 

61 


BETA [13,53 =BETA [11/53 

62 


BETA [7/63 =125. /54. 

63 


SETA [8/63 =-2*/9. 

6 4 


BETA [9/63 --107./9. 

6 5 


3ETA [10/63 *311. /54. 

6 6 


BETA [1 1/63 =-301 ./82. 

A 7 


SETA [12/ 63 =-6./41 . 

68 


BETA [13/ 63 =-289. /82« 

6 9 


'BETA [8/ 73 = 13*/900. 

70 


3ETA [9/ 73 = 67 • /90 • 

71 


3E T A [10/ 73 =-19. /60. 

72 


SETA [1 1/73 =2133. /4100. 

73 


SETA [12/ 73 =-3./205. 

74 


3E t A [13/ 73 =2193. /4100* 

7~ 


BETA [9/ 83 =3. 

76 


BETA [10/ 83 = 17. /6. 

77 


BE t A [1 1 / 83 =45. /82. 

7& 


BETA [12/83 =-3./41. 

7? 


SE t A [13/ 83 =51 */82. 

80 


SETA [1C/93 =-l./12. 

Pi 


SETA [1 1/ 93 =45./164 . 

82 


BETA [12/93 = 3 • / 4 1 . 

8 3 


3ETA [13/ 93 =33./164 . 

84 


BETA [11, 103 =18. /41 . 

p5 


BETA [12/103 =6. /41 . 

86 


BETA [13/ 103 =12. /41 . 

P7 


BETA [13/123 = 1 . 

88 


NS = C 

89 


\R = 0 

/- r> 


N S T = 0 

91 


NRT = C 

52 


PRINT 800 

9 3 

8 C 0 

FQR«AT [1H1///51X/ 14HIMTIAL VALUES! 

54 


CALL PRINT [TI,X/ PPT1 / DTP 1 / DTG I / T0L3 

35 


TsTI 

3 6 


stfp=fpti 

97 


CTS=FPT1-T 

98 


IF r ASS [DTG3 -ABS [DT3I3 3 6/6/7 

99 

7 

DTG=CTGI 

100 

t 

N SF = 0 

in 


NRP = 0 

1C2 


IF [ABSITF-TI3 -ABS[FPT1-TI33 112/121/121 

1 C 3 

11? 

STF p =TF 

104 


DTG=DTG I 
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TABLE A-l. (Continued) 


105 

121 

CALL INTEGRCT/STEP/DTG/TOL/X, 2/KT /NSF/NRF] 

106 


IFtTF-STEP] 161,151,161 

1C7 

161 

T« STEP 

108 


STFP-T+DTP1 

1C9 


I F CABS CDTG] - ABS CDTP1 ] J 8/ 8/ 9 

no 

9 

dtg-dtpi 

111 

8 

IF CABS CTF-T]- ABS [DTP 1]] 132/1*3/1*3 

112 

132 

STFP-TF 

113 


IF CABS CDTG] - ABS CTF/.T3 3 1*3/ 1*3/2024 

11* 

202* 

OTG/TF-T 

115 

1*3 

PRINT 190/NSF/NRF 

116 

190 

FORMAT C 1 H / //// 48X/ 19HINTERMEDI ATE VALUES/3X/ I4/IX/19HG00D 

117 


1AKEN , /1X/I4/1X/15HBAD STEPS TAKEN! 

118 


CALL PRINT [T ,X, FPTl/ DTP1 / DTG I / TOL) 

119 


NS.NS+NSF 

120 


NR»\R+NRF 

1?1 


GO TO 121 

122 

151 

NS«NS+NSF 

123 


NR*NR+NRF 

12* 


NST*NST+NS 

125 


NRT.NRT+NR 

126 


PRINT 840/NST/NRT 

127 

8*0 

FORMAT C1H0/50X, 14HFINAL VALUES *I*/16HG0OD STEPS TAKEN/I*/ 

128 


115HBAD STEPS TAKEN] 

129 


call printctF/x, fpti /Dtpi/Dtgi/TOL] 

130 


GO TO 20 

131 

11 

F0RMATC3E18.il] 

132 


END 


STEPS T 


CeMMgs. ALLOCATION 


77746 ALPH 

77256 BETA 

7722^ CH 


program ALLOCATION 



00011 X 

00015 KT 

00016 I 

00017 J 

00020 NS 

00021 NR 

00022 NST 

00023 NRT 

00024 NSF 

00025 NRF 

00026 T I 

00030 DTOI 

00032 TOL 

00034 TF" 

00036 FPTI 

00040 DTP1 

00042 T 

00044 STEP 

00046 DTG 



SUBPROGRAMS REQUIRED 


print ABS INTEGR 

THE END 
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TABLE A-l. (Continued) 


* 1 SUBROUTINE INTEGRCTI#T,OTS#T0L /X/N/KT/NS/NRJ 

- 2 D I M E N$ I BN XC 2 3/XEC2 3 

« 3 M T * C 

« 4 NS»r 

* 5 NR«C 

• 6 D T 0 = D T S 

• 7 T0=TI 

■ ? 20 xecn»x'ci 3 

« 9 XErE3*yC23 

• 1C STEP-TO+OTG 

* 11 Call RK713CTG/STEP/DTG/T0L vX/ 2 / MS/MR#XE/2 3 

■ 12 TO«STEP 

« 1 3 NS«MS+NS 

= : t NR»NR+MR 

» ;= DTS»OTG 

- 16 NT*NT+MS+MR 

» 1 7 IF CSTEP-T3 240/230/240 

» 18 240 IFCA3SCDTG3 -aBSCT- STEP] 3 210, 210/ 260 

« IS £60 DTG=T-STEP 

■ 20 210 IF CMT-KT3 20/220/220 

« -1 220 T » STEP 

« ?2 230 RETURN 

» 22 END 


PROGRAM ALL0CAT I BN 

DUMMY X 00030 XE 00034 INTEGR 00035 NT 
DUMMY NS DUMMY NR 00036 MS 00037 MR 
DUMMY <T 00040 DTG DUMMY DTS 00042 TO 
DUMMY T I 00044 STEP DUMMY T0L OUMMY T 


subprograms required 

RK71 3 ABS 

THE END 
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TABLE A-l. (Continued) 


• 1 SUBROUTINE RK713 CTI/TF/DT /TOL /X/N/ NS/NR/XE/M3 

■ 2 c seventh order runge-kutta integration with stepsize control 

■ 3 C TF CAN BE GREATER THAN TI OR LESS THAN TI AND RK713 WILL WORK 

• 4 C NS IS THE NUMBER OF SUCCESFULL STEPS TAKEN 

■ 5 C NR IS THE NUMBER OF REJECTED STEPS TAKEN 

■ 6 C N TS THE NU m BER OF DIFFERENTIAL EQUATIONS 

■ 7 C Kt IS MAX number OF iterations 

• 8 C ARRAY F STORES the 13 EVALUATIONS OF THE DIFFERENTIAL EQUATIONS 

• 9 C SUBSCRIPTS FOR ALPHA, BETA/ AND CH ARE +1 GREATER THAN FEHtBERGS 


10 

C 


p[03 IN FEHLBERGS report is IN f ci , j3 


11 

c 


FCIJ IS IN FCI+1,J] 


12 

c 


FEHLBERGS REPORT REFERENCED IS NASA TR R-287 

1 3 

c 


parameters for deq subroutine must be 

STORED IN COMMON 

14 

c 


dimensions must agree with number of 

DIFFERENTIAL EQUATIONS AND 

15 

c 


number of constants in the particular 

FEHLBERG FORMULA USED 

16 



dimension FC13, 2 3/XDUMt 2 3, T EC 2 J 

/ALPHC13J/ 

17 



1BETAC13/123/XC 2 I # CH C 1 33 /XEC2 1 


18 



Common alph,beta,ch 


19 



t*ti 


20 



NS-0 


21 



NR«0 


23 


20 

call deq cx,t,te3 


23 



DO 30 I ■ 1 / N 


24 


30 

F C 1 < II -TECH 


25 



DO 70 K.2,13 


26 



DO 40 I a 1 * N 


27 


40 

xDUMtnaxm 


28 



NN«K-1 


29 



DO 50 I * 1 / N 


30 



DO 50 Jal,NN 


31 


5C 

XDUM C 13 «XDUM CI3 +DT*BETA CK, j3 *F Cj, 13 


32 



TDUM-T+ALPHCK3»DT 


33 



call deq cxdum,tdum,te3 


34 



DO 60 I ■ 1 / N 


35 


60 

FCK/13-TECn 


36 


70 

CONTINUE 


37 



ER«0* 


38 

c 

M 

IS AN INPUT value WHICH DETERMINES the 

NUMBER OF VARIABLES used 1 

39 

c 

THE ERROR CONTROL LOOP 


40 

c 

XE 

IS AN INPUT VECTOR with DIMENSION m WHICH is USED to NORMALIZE 

41 

c 

THE TRUNCATION ERROR COMPUTATIONS IN the 

error control loop 

42 



DO 120 I ■ 1 / M 


43 

140 

TEn3»DT*tFCl,t3+FtU/I3-FC12,T3-FC13, 

, 133*41. /840./XECI3 

44 



If CABSCTECI33-ER) 120/120/130 


45 


130 

ER-ABS CTEI 13 3 


46 


120 

CONTINUE 


47 



DT1*DT 


48 



AK» • 8 


49 



IFCER3 141/142/141 


50 

142 

DT«10.*DT1 


51 



GO TO 150 


52 


10 

OTa CSQRT CSQRT CSORT CTOL/ERI 1 3 3 


53 



0T»AK*DT»DT1 
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TABLE A-l 


(Continued) 


54 


Ir CFR -T0L3 150 a 150a 180 

55 

150 

TF*T+0T1 

56 


00 90 I«1aN 

57 


00 90 U* 1 a 1 3 

58 

90 

X Cl 3 *XCI3+OT1*CHCL3 »FCLaI3 

59 


N S * -N S + 1 

60 


30 T9 230 

61 

180 

NR*NR+ 1 

68 


TF«T 

63 

230 

RETURN 

64 


Eno 


CCMMon ALLBCaT I On 

77746 ALPH 7 7 256 0F.TA 


77224 CH 


PR0GRa m ALL0CAT I 0 k 


0C037 F 

00123 XDUM 

00187 TE 

DUMMY X 

DUMMY XE 

DUMMY NG 

Ouh m y NR 

00133 I 

Dummy [v, 

C0134 K 

00135 NN 

00136 J 

DUM^Y M 

CC137 L 

00140 RK 7l 3 

D01A2 T 

DU-f'Y tI 

DUMMY DT 

00144 TDUM 

00146 ER 

00150 0T1 

C0152 A< 

DUMMY T0L 

DUMMY TF 


SUBPRCGR A-S REQUIRED 


DEQ AbS SORT 

THE Ef-0 


» 1 Su d PT'JTINE PRINT [ T / X / FPT,DTPaDTC-»T0L3 

« 2 dimension xc23 

« 3 print i,t,fpt,dtp,dtg,t0L 

• 4 1 a x 1 1 3 a X [23 

• 5 1 P0P M AT C1H0/5HT «E18.11»2Xa5HFPT «E1 8. 1 1 a 2xa 5HDTP «E18.11a2Xa 

■ 6 15HPTG .E18.11/2X,5HT8L *E 1 8 • 1 1 

• 7 2*/,6H X C13 -ElS. U/2X.5HX C23 »E18. 113 

• 8 RETURN 

• 9 EnD 


PR0GRaM aLLBCaT I 0N 


DUMMY 

X 

00014 

PRINT 

dummy 

T 

DUMMY fpt 

DUMMY 

Twe end 

DTP 

Dummy 

OTG 

dummy 

T0L 
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TABLE A-l. (Continued) 


1 

2 

3 

4 

5 

6 


SUBR0UTINE DEQCX,T,DX3 
DIMENSI9N X [23 > DX [23 

DX C13 *X [13 **2*X [23 
DXC23 »*1 «/X C13 

return 

END 


PR0GRAM ALLOCATION 

DUMMY X DUMMY DX 000Q6 DEO 

THE END 


AASSIQN BI-MT1. 
4REWIND MT1 • 
AF0RT10AD BIU* 



TAB LE A -1 . ( C ontinued) 


NAME 

ENTRY ORIGIN 

LAST 

SIZE/10 

COMMON 

•program 

03507 ?3<*77 

10636 

2656 

17063 


INITIAL VALUES 


T 

•OCOOOCOO°OOE 

00 

FPT ■ 

*50000000000E 10 

DTP • 

.50000000000E 00 

DTG 

• . 10000000000E 00 

TOt • 

•10000000000E oo 

X [1) * 

.10000000003E 

01 

XC2)« 

•10000000000E 01 












final 

VALUES *G90D 

STEPS 

TAKEN OBAD STEPS ' 

TAKEN 


T 

.ACOOCCOOOOOE 

01 

FPT . 

.50000000000E 10 

DTP ■ 

•50000000000E 00 

DTG 

• .10000000000E 00 

TOl • 

tlOOOOOOOOOQE 00 

X[l). 

• 65333167*1 7E 

02 

XC2)> 

*211 9633783 3E -01 









TABLE, A -1. (Continued) 


*&• 

to 


INITIAL VALUES 


T ■ 

•ooooooooopoe 

00 

FPT • 

•50000000000E 

10 

OTP • 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

in 

• 

00 

DTG 

• .10000000000E 00 

T6l ■ 

X[t3. 

.IOOOOOOOOOOE 

01 

Xt23« 

•IOOOOOOOOOOE 

01 













FINAL VALUES 5G080 STEPS 

taken obad steps 

TAKEN 

T • 

f *OOOOOOOOOOE 

01 

FPT • 

.500000000002 

10 

DTP • 

Ui 

o 

o 

o 

a 

o 

o 

o 

o 

o 

o 

in 

• 

00 

DTG 

* .iooooooooooe oo 

T8L • 

Xtl). 

.553682231 22E 

02 

Xt2l» 

•1855593A1 17E' 

-01 














INITIAL VALUES 





T • 

•oooooooooooe 

00 

FPT • 

.50000000000E 

10 

DTP • 

•50000000000E 

00 

DTG 

• t IOOOOOOOOOOE 00 

t»l • 

X [13 • 

•iooooooooooe 

01 

X C23 • 

•IOOOOOOOOOOE 

01 













FINAL 

VALUES 6G08D 

STEPS 

TAKEN OBAD STEPS TAKEN 

T . 

•aooooooooooe 

01 

FPT • 

•50000000000E 

10 

DTP • 

•EOOOOOOOOOOE 

00 

DTG 

« * IOOOOOOOOOOE 00 

tbl ■ 

X{13» 

•5AA73A29571E 

02 

X C23 • 

•18272071378E- 

-01 








•IOOOOOOOOOoEkOI 


.IOOOOOOOOOOE. 01 


.lOOOOOOOOOOE-O* 


•lOOOOOOOOOOE.Ot 



TABLE A-l. (Continued) 


INITIAL VALUES 

T . •OOOOOOOOOCOE 00 FPT ■ »50000000000E 10 DTP • .50000000000E 00 DTG f • 10000000000E 00 TBl • • 10000000000E*03 

X c 13 • .lOOOOCOOOOOE 01 X C23 » . 100000000C0E 01 

FINAL VALUES 8G8BD STEPS TAKEN OBAD STEPS TAKEN 

T « .AOOOOOOOOOOE 01 FPT • •50000000000E 10 DTP • .50000000000E 00 DTG • • 10000000000E 00 T»l ■ t *OOOOOOOODOE*03 

XC1]- .5A5AC130A39E 02 X C2] ■ . 18295966039E-01 

INITIAL VALUES 

T • .OOOOOOOOOOOE 00 FPT • ,50000000000E 10 DTP • .50000000000E 00 DTG f . 10000000000E 00 TOL • • lOOOOOOOOOOE-O* 

X[1J« •lOOOOOOOOCOE 01 XC2]« • 10000000000E 01 

FINAL VALUES 9G08D STEPS TAKEN OBAD 9TEPS TAKEN 

T • .AOOOOOOOCOOE Ol FPT ■ .50000000000E 10 DTP • •50000000000E 00 DTG ■ • 10000000000E 00 TOL • *10000000000E»0* 

X (!) ■ .SASSSSEAATGE 02 XC2)» »183l 1318*38E-01 


00 



& 


TA BLE A -1. ( C ontinued) 


INITIAL VALUES 


T • 

•OOOOOOOOOOOE 

00 

fpt - 

#50000000000E 10 

dtp • 

•50000000000E 

00 

DTQ 

• • 10000000000E 00 

T8l ■ 

•10000000000E»0» 

xc 11 - 

tlOQOOOOOOOCE 

oi 

XC21 ■ 

«10000000000E 01 













FINAL 

VALUES 12G00D 

STEPS 

taken obad steps 

taken 


T • 

• AOOOOOOOOOOE 

Oi 

fpt • 

•50000000000E 10 

DTP • 

.50000000000E 

00 

DTQ 

• • 10000000000E 00 

T8l • 

•10000000000t*05 

xcn- 

•5A59598l?73E 

02 

X t2] ■ 

• U314910017E-01 

w 













INITIAL VALUES 





T 

•OOOOOOOOOOOE 

00 

FPT ■ 

•50000000000E 10 

DTP ■ OOOOOOOOOOOE 00 

DTG 

• .10000000000E 00 

T8L * 

•iooooooooooc»o* 

xcn« 

•10000000000E 

01 

X C23 • 

•10000000000E 01 











FINAL VALUES 15G08D 

STEPS 

TAKEN OBAD STEPS 

TAKEN 


T . 

•AOOOOOOOOOOE 

oi 

fpt • 

•50000000000E 10 

DTP « .50000000000E 00 

DTQ 

T • 10000000000E 00 

T8L • 

•iooooooooooe«o* 

xcn« 

•54597829A17E 

02 

X (23 • 

• 1831 5531 2 42E *01 








TABLE A-l. (Continued) 


INITIAL VALUES 


T 

•OOOOCOOOOOOE 

OC 

FPT • 

•50000000000E 10 

OTP • 

•50000000000E 00 

OTQ 

9 * 10000000000E 00 

Tot • 

#10000000000E-07 

XC»* 

•lOOOOQOOQOQe 

01 

XC21 • 

•10000000000E 01 












PINAL 

VAwUES 2OG08D ! 

STEPS 

taken osad STEPS 

TAKEN 


T 

.^OOOOOOOOOOE 

01 

FPT • 

,50000000000E 10 

OTP « 

►50000000000E 00 

DTG 

• •10000000000E 00 

TQl • 

•10000000000E*07 

X[i}« 

.5*598099<»98E 

02 

X C2] » 

' .18315621928E-01 














INITIAL VALUES 






T 

•oooooooooooe 

00 

FPT • 

•50000000000E 

10 

DTP • 

•sooooooooooe 

00 

DTG 

* ^ tlOOOOOOOOOOE 00 

T&U • 

.lOOOOOOOOOQE-Ot 

XtD* 

.lOOOOCOOOOOE 

01 

X [?) * 

tlOOOOOOOOOOE 

01 














final 

VALUES 26G60D ! 

STEPS 

TAKEN 08AD STEPS TAKEN 


T . 

.AOOOOOOOOOOE 

01 

FPT ■ 

•50000000000E 

10 

DTP • 

•50000000000E 

00 

DTG 

* tlOOOOOOOOOOE 00 

T8L * 

•lOOOOOQOOOoE.Ot 

xcn- 

•5A598133849E 

02 

X 12) ■ 

• 1831 5633457E 1 

•01 










TABLE A-l. (Continued) 


INITIAL VALUES 


T 

.OOOOOOOOOOOE 00 

FPT . 

•50000000000E 10 

DTP • 

•50000000000E 00 

DTQ 

• tlOOOOOOOOOOE 00 

T0L • 

• lOOOOOOOOOOExO. 

v xt»« 

.looooooooooe oi * 

XC23- 

-•10000000000E 01 











FINAL 

VALUES 30GB0D 1 

STEPS 

TAKEN OBAD STEPS ' 

TAKEN 


T 

.OOOOOOOOOOOE 01 

FPT • 

,50000000000E 10 

DTP * 

.SOOOOOOOOOOE oo 

OTG 

- tlOOOOOOOOOOE 00 

TOL • 

.lOOOOOOOOOOE.O* 

XC»« 

• 50598 135326C 02 

XC23- 

■ *18315633952E-01 








INITIAL VALUES 

T . •OOOOOOOOOOOE 00 FPT • . 50000000000E 10 DTP • .50000000000E 00 DTQ • . 10000000000E 00 T&L • *10000000000E*lO 

K£l3 • .10000000000E 01 K C23 • • 10000000000E 01 

FINAL VALUES MG08D STEPS TAKEN OBAD 9TEPS TAKEN 

T • .OOOOOOOOOOOE 01 FPT i .SOOOOOOOOOOE 10 OTP • •SOOOQOOOOOOE 00 DTQ • * 10000000000E 00 TOl » *10000000000E«il 

*£!)• .5A59B13300AE 02 XC2J- • IS3194331 70C*0l * 

; S ' 


TABLE A-l. (Continued) 


INITIAL VALUES 


T 

•OOOOOOOOOOOE 

00 

FPT • 

tSOOOOOOOOOOE 

00 

OTP • tSOOOOOOOOOOE 

00 

DTG 

• .10000000000E 

00 

T0L p tlOOOOOOOOOOE-lO 

X[l]« 

tlOOOOOOOOOOE 

01 

X C23 • 

tlOOOOOOOOOOE 

01 












INTERMEDIATE VALUES 

6 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

tSOOOOOOOOOOE 

00 

FPT « 

•50000000000E 

00 

DTP ■ .50000000000E 

00 

DTG 

• .10000000000E 

00 

TOL ■ *10000000000E*10 

xcn« 

•16A8721270AE 

01 

XC23- 

•60653065960E 

00 












INTERMEDIATE values 

7 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

tlOOOOOOOOOOE 

01 

FPT • 

.50000000000E 

00 

DTP « .50000000000E 

00 

DTG 

■ . 10000000000E 

00 

TOL ■ .10000000000E-10 

xcn- 

•27182818P65E 

01 

XC2] • 

•367879AA091E 

00 












intermediate values 

10 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T • 

.15000000000E 

01 

FPT ■ 

.50000000000E 00 DTP • tSOOOOOOOOOOE 

00 

DTG 

• 1 10000000000E 

00 

TOL • .lOOOOOOOOOOE.iO 

xcn. 

•A4816890A00E 

01 

XC23- 

.22313015963E 00 










INTERMEDIATE VALUES 

9 

GOOD 

STEPS TAKEN t 

0 

BAD STEPS TAKEN 

T 

tSOOOOOOOOOOE 

01 

FPT • 

tSOOOOOOOOOOE 00 DTP ■ tSOOOOOOOOOOE 

00 

DTG 

• tlOOOOOOOOOOE 

00 

T0L • tlOOOOOOOOOOE- 10 

xcn. 

•73890560SC5E 

01 

XC23 ■ 

tl353352823AE 00 










intermediate values 

8 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

•25000000000E 

01 

FPT • 

.50000000000E 00 DTP - .50000000000E 

00 

DTG 

* . 10000000000E 

00 

TOl ■ .lOOOOOOOOOOE-10 

xcn- 

t 1 21 82A937A0E 

02 

X C23 ■ 

•8208A997121E-01 










intermediate values 

7 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

.30000000000E 

01 

cPT . 

.50000000000E 00 DTP • .50000000000E 

00 

DTG 

■ tlOOOOOOOOOOE 

00 

TOl ■ tlOOOOOOOOOOE*iO 

xcn. 

.2008S535926E 

02 

X C 2 3 » 

•49787065B82E-01 










intermediate values 

6 

G00D 

STEPS TAKEN , 

0; 

BAD STEPS TAKEN 

T • 

•35000000000E 

01 

FPT ■ 

tSOOOOOOOOOOE 00 DTP ■ tSOlOOOOOOOOE 

00 

DTG 

• tlOOOOOOOOOOE 

00 

TOL ■ tlOOOOOOOOOOE. 10 

xcn- 

.33115A47A75E 

02 

X C23 ■ 

•30197379324E-01 







I 



TABLEA-1 


(Concluded) 


FPT ■ 
X C2J • 


FINAL VALUES 59G90D STEPS TAKEN OBAD STEPS TAKEN 

.50000000000E 00 DTP • •50000000000E 00 DTO • 10000000000E 00 TBl * • lOOOOOOOOOQE»tO 

•18315632132E*01 ' ' 


AOOOOOOOOOOE 01 
5A59S129917E 02 



TABLE A -2. TEST PROBLEM FI 


^ASSIGN S«MT0/SI CR/B8 MT1/L0 LP. 

^REWIND MT1 * 

AF0RTRAN B0/L0* 

1 DIMENSION XC43 

2 DIMENSION ALPH[i33,BETACl3,i23/CH[133 

3 C0MM0N GMM 

4 C0MM0N A|_PH/BETA,CW 

5 20 READ 11, TI/DTGI 

6 READ 11, [XCI3/I-1/43 

7 READ 11,TF/FPT1#DTP1 

8 READ 11/GMM/T0L 

9 READ 2C25, KT 


SUBROUTINE 


10 

2025 

FORMAT C X 43 

11 

C 

CONSTANTS F0R INTEGRAT 

12 


D0 60 1=1/13 

13 


D0 50 J* 1 , 12 

14 

50 

be t a C I , j3 »o. 

15 


alph cn =o* 

16 

60 

CHCI3 -0. 

17 


CH C63 »34./105. 

18 


CH [73 *9./35. 

19 


CH [83 » CH [ 73 

20 


CH [93 = 9 . /280 • 

21 


CH [ 1 03 =CH [93 

22 


CH [ 1 23 =41 ./840. 

23 


CH [ 1 33 =CH [123 

24 


ALPH [23 = 2 • /27 • 

25 


ALPH [33 -1 ./9* 

26 


A|_PH [43 = 1 ,/6» 

27 


ALPH [53 = 5 . / 1 2 • 

28 


AL p H [63 • ,5 

29 


AL D H [73 «5./6* 

30 


ALPH [83 =1 ,/6« 

31 


ALPH [93 » 2 • /3 • 

32 


ALPH [103 ■ 1 ,/3» 

33 


ALPH [113*1* 

3^ 


ALPH [133 »1 . 

35 


BETA [2/ 13 =2./27. 

36 


BETA [3/ 13 « 1 ./36. 

37 


BETA [4/ 13 =l./24. 

38 


BETA [5/ 13 »5./12. 

39 


BETA [6/13 * • O 5 

40 


BETA [7/ 13 --25./108. 

41 


BETA [8/ 13 =31 ./300. 

42 


BETA [9/ 13 =2. 

43 


BETA [10/ 13 =-91 ./108* 

44 


BETA [1 1, 13 = 2383, /4 100. 

45 


BETA [12/ 13 * 3 • /205 . 

46 


BETA C13/ 13 --1777./4100 

47 


BETA C3/ 23 ■ 1 • / 1 2 • 

48 


BETA [4/33 «1 ,/8. 

49 


BETA C5/33 --25. /16. 

50 


BETAC5/43--BETAC5/33 



51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 


TABLE A -2. (Continued) 

BETA C6/ 43 ■ *25 
BETA C 7# 4] * 125. /108,« 

BETA [9/ 43 *-53*/6. 

BETAC10, 43*23. /108. 

BETAC11/43--341./164. 

BETAC13/43-BETAC11/43 

BETA [6# 53 * *2 

BETA C7/ 53 *-65»/27» 

BETA [8/ 53 -61./225. 

BETA [9/53 *704. /45. 

BETA CIO/53 "-976./135. 

BETA [11/53 *4496. /1025. 

BETA C 13# 53 *BETA[ll/53 
BETA C7/63 -125./54. 

BETA [8/63 «-2»/9. 

BETA [9/63 ■-107./9. 

BETA [10/63 *311. /54» 

BETA [11/63 --301./82. 

BETA [12/63 *“6./41* 

BETA[i3/63«-289./82* 

BETA[8/73 -13./900. 

BETA [9/73 *67. /90. 

BETA [10/73 --19./60. 

BETA[11, 73*2133. /4100. 

BETA [12/ 73 *»3./ 205. 
BETA[13 / 73=2193./4100. 

BETA [9/83 *3. 

BETA [10/83 -17. /6. 

BETA [11/83 -45. /82. 

BETA [12/ 83 *-3./41 • 

BETA[1 3/83*51./ 82. 

BETA [10/93 --1./12. 

BETA[ll/93«45./164. 

BETA [12/93 • 3 * /A 1 • 

BETA [13/93 -33./16A, 

BETA[U/103 -18./41* 

BETA [12/103 -6./41. 

BETA[13/103 -12./41, 

BETA [13/ 123 *1. 

NS-0 
NR«0 
NST«0 
NRT*0 
PRINT 800 

800 FORMAT [1H1///51X/14HINITIAL VALUES! 

call print cti,X/ fpti /DTpi/DtgI/TOlj 
t-t I 

step-fpti 

DTG-FPT1-T 

IF [ABS [DTG3 -ABS CDTGI3 3 6/6/7 
7 DTG-DTGI 

6 NSF-0 

NRF*0 

IF CABS CTF-TI3 -ABS CFPTl*TU)il2* 121*121 



TABLE A -2. (Continued) 


105 

106 

112 

107 

108 

121 

109 

no 

in 

161 

112 

9 

113 

3 

11* 

115 

132 

116 

2024 

117 

143 

118 

119 

120 
121 
122 
123 

190 

12* 

125 

126 

127 

128 

151 

129 

130 

131 

132 

840 

133 

13^ 

11 


STEP-TF 

DTG-DTGI ' 

CALL INTEGRCT/STEP/DTG/TBL/X, 4/KT /NSF/NRF! 

IFCTF-STEPH61, 151/161 

T.STEP 

STEP*T+DTP1 

IF TABS CDTG! -ABS CDTP1!! 8/8/9 
DTG-DTPl 

IF CABS CTF-T3 -ABS [DTP 131 132/ 143, 143 
STFP-TF 

IF t ABS CDTG! - ABS CTF-T1! 1 43/ 143/ 202 A 
DTG.TF-T 

PRINT 190, NSF/NRF 

F0RMATC1H ///// A8X/19HINTFRMEDI ate VALUES/ 3X, I A/ 1 X/ 19HG0OO STEPS T 
1AKFN / / 1 X/ 1 4/ 1 X/ 1 5HBAD STEPS TAKEN! 

CALL PRINT tT , X/ FPT1/DTP1/DTGI/T8L) 

NS»NS+NSF 

NR»NR+NRF 

GO TO 121 

NS*NS+NSF 

NR.NR+NRF 

NST»NST+NS 

NRT.NRT+NR 

PRINT 840/NST/NRT 

FORMAT C1H0/50X, HHFINAL VALUES , I4/16HG00D STEPS TAKEN/ 14 , 

115MBAD STEPS TAKEN! 


call printctF/X, 
GO TO 20 
FORMAT C3E1 8 » 1 1 ! 
END 


PPT1 /DTP1/DTGI/T8L1 


COMMON ALLOCATION 

77776 GMM 77744 ALPH 77254 BETA 77222 CH 


PROGRAM ALLOCATION 


00011 X 

00021 I 

00022 KT 

00023 J 

00024 NO 

00025 NR 

00026 NST 

00027 NRT 

00030 NSF 

00031 NRF 

00032 T I 

00034 DTQI 

00036 TF 

00040 FPT 1 

00042 DTP1 

00044 T0L 

00046 T 

00050 STEP 

00052 DTG 



SUBPROGRAMS REQUIRED 

PRINT ABS INTEGR 

the end 


51 



TABLE A-2. (Continued) 


1 


SUBROUTINE INTEflRCTI/T/DTS/TOL /X/ 

2 


DIMENSION Xt ‘4 )*XCU 3 

3 


NT«0 

4 


NS-0 

5 


NR-0 

6 


DTG'DTS 

7 


TO-TI 

8 

20 

R2»XC13**24Xt23**2 

9 


V2»XC33#*2+XC43#*2 

10 


R-SQRTCR23 

11 


V-SQRT CV23 

12 


DO 1 1-1/2 

13 


XECI3-R 

14 

1 

XE 1 1 +2] »V 

15 


STEP-TO+DTG 

16 


Call RK713CT0/STEP/DTG/T8L /X/ 4 / 

17 


TO-STEP 

18 


NS»MS+NS 

19 


nr-nr+mr 

20 


DTS-DTG 

21 


nt»nt+ms+mr 

2? 


IF CSTEP-T3 240/230/240 

23 

240 

IF CABSCDTG3 -ABS IT- STEP3 3 210/230/ 260 

?4 

260 

DTG-T-STEP 

25 

210 

IFtNT-KTJ 20/220/220 

26 

220 

T-STEP 

27 

230 

RETURN 

28 


END 


/X/N/KT/N*,N*J 


MS/MR/XE/4 3 


PROGRAM ALLOCATION 


DUMMY X 

00031 XE 

00041 INTEGR 

00042 NT 

DUMMY NS 

DUMMY NR 

00043 I 

00044 MS 

00045 MR 

DUMMY KT 

00046 DTG 

DUMMY DTS 

0005° TO 

DUMMY T I 

00052 R2 

00054 V2 

00056 R 

00060 V 

00062 STEP 

DUMMY TBl. 

DUMMY t 





subprograms required 

SORT RK713 ABS 

THE end 
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TABLE A-2. (Continued) 


• 1 subroutine RK713 cti,tf,dt ,tol ,x,n, ns,nr,xe,m3 

• 2 c seventh order runge-kutta integration with stepsize control. 

• 3 C TF CAN BE GREATFR THAN Tl OR LESS THAN TI AND RK713 WILL PORK 

• 4 C NS IS THE NUHRrR OF SUCCESFULL STEPS TAKEN 

» 5 C NR IS THE NUHBfR OF REJECTED STEPS TAKEN 

• 6 C N IS THE NUMBER OF DIFFERENTIAL EQUATIONS 

• 7 C KT IS MAX NUMBER OF ITERATIONS 

■ 8 C ARRAY F STORES THE 13 (■ VALUATIONS OF THE DIFFERENTIAL EQUATIONS 

■ 9 C SUBSCRIPTS FOR ALPHA* BETA* AND CH ARE +1 GREATER THAN FEHLBERGS 

• 10 C F [B] JN FEHLBERGS REPORT IS IN Ftl*J3 

• 11 C F[II IS IN FU + hjJ 

• 12 C FEHLBERGS REPORT REFERENCED IS NASA TR R-287 

• i3 c parameters for deq subroutine must be stored in common 

• i4 c dimensions must agree with number of differential equations and 

• is c number of constants in the particular fehlberg formula used 

• 16 DIMENSION F C 13* 4 J*XDUMC 4 ] * TE C 4 ] * ALPH C133 * 

■ 1 7 16ETAC13,121,XC 4 I , CH C 1 31 * XE C4 3 

• 18 Cl) “MON GMM 

• 1 9 CH"MON ALPH, BETA, CH 

• ?0 T.TI 

• 2 1 N5«0 

• 22 NR ■ 0 

• ?3 20 CALL DEQ CX,T,TE3 

• 24 DO 30 I ■ 1 * N 

« ?5 30 F c 1 , I3»TECI3 

« 26 DO 70 K ■ 2 * 13 

• 27 DO 40 I - 1 * N 

• 28 40 XDIJM C 1 3 *X C 13 

■ 29 NN«K-1 

• 30 DO 50 I » 1 * N 

■ 31 DO 50 J»1*NN 

■ 32 50 XD'JM CI3 «XDUMCI3 +DT«9ETA CK* j3 »F Cj* 13 

■ 33 Tdi;m»T + ALPhCK3 * DT 

• 34 call deq cxdum,tdum*te3 

■ 35 DO 60 I ■ 1 * N 

- 36 60 F [K* I] «TE CI3 

• 37 70 CONTINUE 

■ 38 ER ■ 0 • 

« 39 C M IS AN INPUT VALUE WHICH DETERMINES THE NUMBER Or VARIABLES USED IN 

■ 40 c the eR r or CONTROL LOOP 

« 41 C XE is AN Input VECTOR with DIMENSION m which is used to NORMALIZE 

• 42 C THE TRJNCATION ERROR COMPUTATIONS IN THE ERROR CONTROL loop 

« 43 DO 120 I • 1 , M 

• 44 140 TECI3 -DT» CF Cl, 13 ♦FCll* I3-FC12, I3-FC13* 133*41, /840./XECI3 

■ 45 IF CABS CTE C 1 3 3 «ER3 120, 120,130 

• 46 130 ER-ABS CTE CI3 3 

• 47 120 CONTINUE 

■ 48 DTl'DT 

■ 49 A< ■ , 8 

• 50 IF CER3 141, 142, 141 

■ 51 142 DT* 10, *DT1 

■ 52 GO TO 150 

■ 53 141 DT-CSQRTCSQRTCSQRTCTOL/ER3333 
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TABLE A -2. (Continued) 


■ 54 


DT»AK*DT*DTt 




• 55 


IF CER •T8LJ 

180# 180# ISO 



• 56 

150 

TF-T+DT1 




■ 57 


D0 90 I - 1 # N 




■ 58 


DB 90 U-l/13 




■ 59 

90 

X C 1 3 •XCI)+DTl*CHCLJ*FCL*n 



» 60 


NS»NS+1 




■ 61 


GO TO 230 




■ 62 

180 

NR.NR+1 




■ 63 


TF-T 




■ 64 

230 

RETURN 




• 65 


END 




COMMON ALLOCATION 




77776 

GMM 

77744 ALPH 

77254 BETA 

77222 

CH 

program 

allocation 




00037 

F 

00207 XDUM 

00217 TE 

DUMMY 

X 

DUMMY 

XE 

DUMMY NS 

dummy nr 

00227 

I 

DUMMY 

N 

00230 K 

00231 NN 

00232 

J 

DUMMY 

M 

00233 L 

00234 RK713 

00236 

T 

DUMMY 

T I 

DUMMY dt 

00240 TOUM 

00242 

ER 

00244 

DTI 

00246 AK 

dummy TOL 

DUMMY 

TF 


SUBPROGRAMS REQUIRED 

deq abs sort 

the. end 


• 1 SUBROUTINE PRINT CTiXi FPT,DTP/DTG*TOLI 

■ 2 D I HENS I BN XC4J 

» 3 COMMON GMM 

■ 4 PRINT l,T/FPT/DTP,DTG,TOL 

■ 5 1/ CXCn»I«l/4J>GMM 

• 6 1 FORMAT C1H0/5HT -E 1 8 . 1 1 , 2X/ 5HFPT «E1 8. 1 1 * 2X, 5HDTP •E18.H*2X* 

■ 7 15HDTG »E18.11»2X,5HT0L -ElS*!! 

■ 8 2 >/, 6 H Xtll«E18.n#2X/5HXC2J"E18.U#2X#5HXC3J»El8.il*2X#5HXtO»E18i 

• 9 3U/2X,5HGMM *E18*11J 

• 10 return 

• n end 


COMMON ALLOCATION 
77776 GMM 

PROGRAM ALLOCATION 

DUMMY X 00014 I 00015 PRINT DUMMY T 

DUMMY FPT DUMMY DTP DUMMY DT* i DUMMY TH 



1 

2 

3 

4 

5 

6 

7 

8 
9 


10 

11 

12 

13 

14 

15 

16 

17 

18 


TABLE A -2. (Continued) 

subroutine deqcx,t#dx3 

DIMENSION X [43 i DXC43 

C0MM0N GMM 
Cl*GMM-l . 

C2“ X [ 1 3 +GMM 
C3*C2-1 • 

R1?»C2**2+X [23##2 
R22»C3#-*2+X [23 #*2 
R1«SQRT [R123 
R2»SQRT[R223 
DEN1«C1/R1/R12 
DEN2-GMM/R2/R22 
DX [13 «XC33 
DX C23 «XC43 

DXC33»XC13+2t#X [43 +DEN1*C2-DEN2*C3 
DX C43 “X [23 -2*«X[33-k[DENl*DEN23*Xt23 
RETURN 
END 


C0MM0N ALLOCATION 
77776 GMM 

program allocation 

DUMMY X DUMMY DX 00013 DEQ 00015 Cl 

00017 C2 00021 C3 00023 R12 00025 R22 

00027 R 1 00031 R2 00033 DENI 00035 DEN2 


SUBPROGRAMS REQUIRED 

SORT 
THE END 


^ASSIGN BI-MJ1. 
aREWIsD MT1 • 
aFoRTL&AD B I u • 



NAME 

•program 


T • 

Xtl)« * 

T . 
XU) * 


TAB LE A -2 . ( C ontinued) 

entry origin last size/jo common * Base » 

03507 03*77 11237 29)3 17061 


INITIAL VALUES 


•OOOOOOOOOOOE 00 FPT > • 10000000000E 11 OTP ■ .10000000000E 01 DTG • • 10000000000E 00 TOL f 

•99400000000E 00 X C23 • r tOOOOOOOOOOOE' 00 X C33 • *OOOOOOOOOOOE 00 Xt*1- -t203l7326296t 01 QMM • 

FINAL VALUES 160000 STEPS TAKEN 12BAD STEPS TAKEN 


•1112*3*0337E 0? 
# 1*63*41 4560E 0? 


FPT • .10000000000E il 

X C2) • .tl0**795776lE 03 


OTP ' .10000000000E 01 DTG • « 10000000000E 00 TOt f 
XC3)« ••10353360093E 03 Xt*)f •. 2603*09601 IE 02 QMM ' 


•10000000000E 00 
•12277471000E-01 


«10000000000E 00 
•1 2277*71 OOOE«Ol 


TABLE A -2. (Continued) 


INITIAL VALUES 


T . 

♦OOOOOOOOOOOE 

00 

FPT • 

•10000000000E 

11 

OTP • 

.ICOOOOOOOOOE 

01 

DTG • 

•10000000000E oo 

T8l * 

•10000000000E«01 

*tn. 

•99AOOOOOOOOE 

00 

x cS] ■ 

•OOOOOOOOOOOE 

00 

XC33- 

•ocoooooooooE 

00 

XC*Jf 

-•20317326296E 01 

GMM ■ 

•12277A71000E-01 







EINAL 

VALUES 21G88D : 

STEPS TAKEN 9BA0 STEPS ' 

TAKEN 


T • 

•1 1 12*3*0337E 

02 

FPT . 

•10000000000E 

11 

DTP ■ 

Ui 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

01 

DT3 • 

ilOOOOOOOOOOE 00 

T8L ■ 

.lOOOOOOOOOOE-Ol 


.A609AQ20?03e 

00 

X [23 • 

*t52836!97765E 

o 

o 

XC3). 

•6T6S6511117E 

00 

X E*J • 

-•36212622822E-01 

GMM * 

«12277*7l000E*0l 


I 

5 

■ | 

■ 

initial values = 


T 

.OOOOOOOOOOOE 

00 

FPT • 

,10000000000E 11 

DTP ■ 

•ICOOOOOOOOOE 01 

OTG • 

ilOOOOOOOOOOE 

00 

T8L • 

.iooooooooooE-oe 

xcn- 

•99AOOOOOOOOE 

00 

X C2J » 

•OOOOOOOOOOOE oo 

X [33 ■ 

•OOOOOOOOOOOE 00 

xc*)» 

'••20317326296E 

01 

GMM ■ 

•12277871000E.01 






FINAL 

VALUES 380980 ! 

STEPS TAKEN 18BAD STEPS ' 

TAKEN 


T . 

.1112*3*03372 

02 

FPT • 

•10000000000E 11 

DTP « 

•10000000000E 01 

OTG ■ 

IlOOOOOOOOOOE 

00 

T8L - 

•10000000000E-02 

X(l)» 

•9939A287931E 

00 

Xt2)» 

•206750939S9E-03 

Xt31" 

•30556134955E-01 

X [A] ■ 

•t20396719553E 

01 

GMM • 

•12277471000E.01 




cn 

<1 



TABLE A -2. (Continued) 


T . 

.OOOOOOOOOOOE 

00 

FPT - 

, lOOOOOOOOOOE 

11 

X Cl] - 

•99400000000E 

00 

X C2] • 

•OOOOOOOOOOOE 

00 

T • 

• U124340337E 

02 

FPT • 

.lOOOOOOOOOOE 

11 

1(1]* 

.99395953273E 

00 

XC23- 

.31208898899E- 

-03 


T . 

•OOOOOOOOOOOE 

00 

FPT • 

.lOOOOOOOOOOE 

11 

• Ill • 

•99400000000E 

00 

* X C21 • ‘ 

’ .OOOOOOOOOOOE' 

00 

T . 

• U12A340337E 

02 

FPT i 

.lOOOOOOOOOOE 

11 

Xtll* 

•99398586886E 

00 

X C23 • 

• .5272500431 oE* 

•04 


INITIAL VALUES 

OTP • .lOOOOOOOOOOE 01 DTG • . 10000000000E 00 T8L p 

X £33 • .OOOOOOOOOOOE 00 X * -• 2Q317326296E 01 SMH » 

FINAL VALUES 49G88D STEPS TAKEN 23BAD STEPS TAKEN 

OTP • .lOOOOOOOOOOE 01 DTG * . lOOOOOOOOOOE 00 T8L i 

Xt33- .47930058035E-01 X ■ *. 2Q362226800E 01 GUM p 


INITIAL VALUES 

DTP • .lOOOOOOOOOOE 01 DTG ■ # lOOOOOOOOOOE 00 T8u ■ 

X C33 • .OOOOOOOOOOOE 00 * X tM • 1 -.20317326296E 01 fiMM ■ 

FINAL VALUES 59G88D STEPS TAKEN 260AD STEPS TAKEN 

DTP ■ .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 T8L p 

X C 3 3 • * .85866 7946 34E *02 X [4] T ' •* 2Q338573144C 01 GHM f 


* lOOOOOOOOOOE -03 
•12277471000E-01 


•10000000000E*03 
. 1 227747 1000E*01 


•iooooooooooe-o^ 

•12277471000E-01 


.1OOOOOO0OOOE«O4 

•1227747IO00E-OI 



TABLE A-2. (Continued) 


INITIAL 'VALUES 


T 

.OOOOOOOOOOOE 

00 

FPT . 

.lOOOOOOOOOOE 11 

DTP • 

•lOOOOOOOOOOE 01 

DTG 

■ • lOOOOOOOOOOE 00 

tbl ■ 

•1OOOOOOOOO0E-O5 

X [13 ■ 

.99400000000E 

00 

XC23- 

•OOOOOOOOOOOE 00 

XC33- 

•ooooooooooot 00 

X(A3 

■ -•20317326296E01 

GMM * 

•12277*710002-01 






F INAL 

VALUES 71G80D ! 

STEPS 

TAKEN 29BAD STEPS 1 

TAKEN 


T . 

•1112A3A0337E 

02 

FPT . 

•lOOOOOOOOOOE 11 

DTP • 

•lOOOOOOOOOOE 01 

DTG 

• • lOOOOOOOOOOE 00 

tol ■ 

•IOOOOOOOOOOE-05 

xcn« 

.99399278977E 

00 

X C23 » 

-•231A3635052E-0A 

X C33 * 

-.379643A0325E-02 

X[4) 

? -.2Q328307961E 01 

GMM * 

• 1227747l000E«0l 


INITIAL VALUES 

T • *OOOOOOOOOOOE 00 FPT ■ . 10000000000E 11 DTP • . 10000000000E 01 DTG • t lOOOOOOOOOOE 00 T8L • • 1000000000QE-06 

Xtl)« .99AOOOOOOOOE 00 XC23- ’ .OOOOOOOOOOOE 00 X C33 - .OdOOOOOOOOOE 00' X [*) • •• 20317326296E 01 GMM * t i2277471000E*Ol 

FINAL VALUES 88GOOD STEPS TAKEN 288AD STEPS TAKEN 

T ■ •111243%o337E 02 FPT ■ .lOOOOOOOOOOE tl DTP • .lOOOOOOOOOOE 01 DTG • » 10000000000E 00 T»L ■ • 10000000000E*06 

XtlU .99399965197E 00 XC23« - . 18169*921 85E-05 X C33 - - .29u5260859E*03 X [43 •“ -t203l786166lE 01 GMM • *1227747*OO0E-Ol 


C7I 

CD 



TAB LE A -2 • ( C ontinued) 


T 

•OOOOOOOOOOOE 

00 

FPT • 

•lOOOOOOOOOOE 

11 

xcn • 

•99400000000E 

00 

XC21 ■ 

•OOOOOOOOOOOE 

00 

T > 

•1112A3A0337E 

02 

FPT • 

•lOOOOOOOOOOE 

11 

Xtl3« 

•99399987755E 

00 

X C23 ■ 

-•5286315A165E- 

•06 


INITIAL VALUES 

DTP - .10000000000E 01 DTG ■ . 10000000000E 00 TOL • 

X C33 ■ .OOOOOOOOOOOE 00 *XCA)> -•20317326296k Oi GMM • 

FINAL VALUES 110G00D STEPS TAKEN 23BAD STEPS TAKEN 

OTP ■ •lOOOOOOOOOOE 01 DTG • * 10000000000E 00 TOL • 

X [33 ■ -.85A80979A0SE-0A X C43 f - 1 2031 751 M50E 01 GMM » 


INITIAL VALUES 

T > •OOOOOOOOOOOE 00 FPT ■ , 10000000000E 11 DTP • • 1 OQOOOOOOOOE 01 DTG • • 10000000000E 00 TOL ■ 

XtUm .99400000000E 00 XC2}«' •OOOOOOOOOOOC 00 * X C33 • " *00600000000E 00 ' X W ■ 20317326296E 01 GMM • 

FINAL VALUES 1A2G00D STEPS TAKEN 17BAD STEPS TAKEN 

T • •U12A3A0337E 02 FPT ■ , lOOOOOOOOOOE U DTP ■ .lOOOOOOOOOOE 01 DTG ■ •lOOOOOOOOOOE 00 TOL • 

XCl)* .99399992998E 00 X[2]« - * 25967018380E-06 XC3]t • . *2320S028()5E-0A Xt*J« ••20317A335A6E 01 GMM ■ 


•10000000000E-07 

•1227747l000E-0t 


•10000000000E-07 

•12277471000E*01 


•10000000000E*OS 
• 1227747l000E*0l 


*10000000000E«N 
• 12277A7I000E*0I 


TABLE A-2. (Continued) 


initial values 


T 

•0QQ0000000QE 

00 

FPT « 

•lOOQOOOOOOOE 

11 

DTP • 

• 10000000000E 

01 

DTG - 

•10000000000E 00 

TOL 

■ 

XlU* 

•99400000000E 

00 

XC23i 

•oooooooooooe 

00 

X C3) * 

•OOOOOOOOOOOE 

00 

XCAl- 

-♦20317326296E 01 

gmm 

■ 







FINAL 

values 187GO0D « 

STEPS TAKEN HBAD STEPS TAKEN 


T • 

•1112A340337E 

02 

FPT - 

•iooooooooooe 

11 

DTP ■ 

•10000000000E 

01 

DTG ■ 

.10000000000E 00 

T0L 

■ 

xc». 

•99399990570E 

00 

X C2] ■ 

••3385778 82 84E- 

■06 

X C31 ■ 

• *55281 726A1 8E' 

•0* 


■•20317A70697E 01 

gun 

• 







INITIAL VALUES 





T 

•oooooooooooe 

00 

FPT ■ 

•lOOQOOOOOOOE 11 

DTP ■ 

•lOOOOOOOOOOE 01 

DTG • 

•10000000000E 00 

T0L 

■ 

XClJ- 

•99A00000000E 

00 

XT23- 

•OOOOOOOOOOOE 'DO 

X [33 ■ 

•OCOOOOOOOOOE 00 

‘xr*i ■ 

••20317326296E 01 

GMH 

• 






FINAL 

VALUES 284G00D ! 

STEPS TAKEN 6BAD STEPS TAKEN 


T . 

• U12A3A0337E 

02 

FPT . 

.10000000000E 11 

DTP • 

•10000000000E 01 

DTG ■ 

•lOOQOOOOOOOE 00 

T0L 

• 

xtn- 

•99399975A5BE 

00 

X C23 ■ 

-.89878A53735E-06 

X t3) » 

-.U661985150E-03 

XC4l« 

**2031770196 9E 01 

gmm 

• 


•IOOOOOOOOO0E-O9 

tl227747!000E-0l 


•lOOOOOOOOOOE-Of 
• 1227747100oE»Ol 


• iooooooooooemo 
•12277471000E-01 


•10000000000E*10 
• l227747lOQOE»Ol 



TABLE A -2. (Continued) 

INITIAL VALUES 


t • .oooooooooooe 

00 

FPT ■ .lOOOOOOOOOOE 

01 

DTP • 

•lOOOOOQOOOOE 

01 

DT3 • .lOOOOOOOOOOE 

00 

T8L r 

• lOOOOOOOOOOE- io 

K CU • «994oc0000.00£ 

00 

X C2} • .OOOOOOOOOOOE' 00 

X 1 33 • 

•OOOOOOGOOOdt 

00 

X C*3 ■ ••20317326296E 

01 

gmm • 

. 1227747lOOOE*Ol 




INTERMEDIATE VALUES 

91 

GOOD STEPS TAKEN $ 

3 l 

BAO STEPS TAKEN 

T - ♦ lOOOOOOOOOOE 

01 

FPT ■ .lOOOOOOOOOOE 

01 

DTP • 

•lOOOOOOOOOOE 

01 

DTG ■ .lOOOOOOOOOOE 

00 

TOL • 

.lOOOOOOOOOOE. 10 

ttUu .51306969334E* 

01 

XC23- .33261806494E 

00 

X C3] • 

• . 1 7725293980E 

01 

X [A] ? .30780989933E 

00 

GMM • 

« 1227747l000E*0l 




INTERMEDIATE VALUES 

25 

GOOD STEPS TAKEN $ 

0 BAO STEPS TAKEN 

T • •20000000000E 

01 

FPT • .lOOOOOOOOOOE 

01 

DTP ■ 

•10000000000E 

01 

DTG ■ .lOOOOOOOOOOE 

00 

TOL • 

.10000000000E-10 

X (13 * **93l3i060560E 
• 

00 

XC21- .29597656687E 

00 

XC3J. 

••9312915S899E 

00 

XCAJf .31006193716E 

00 

GMM • 

.12277A71000E-0* 




INTERMEDIATE VALUES 

10 

GOOD STEPS TAKEN * 

0 

BAO STEPS TAKEN 

t • .30000000000E 

01 

FPT ■ .lOOOOOOOOOOE 

01 

DTP ' 

•10000000000E 

01 

DTG ? .lOOOOOOOOOOE 

00 

TGL * 

. lOOOOOOOOOOE* 1. 

X(l]» ••9761#41234%E 

00 

Xt23» .77776335199E 

oo 

X (33 ■ 

•33149953230E 

00 

XC*3f .565A0311203E 

00 

GMM ■ 

.12277471000E-.S 




INTERMEDIATE VALUES 

8 

GOOD STEPS TAKEN * 

0 

BAO STEPS TAKEN 

T • .40000000000E 

01 

FPT • .lOOOOOOOOOOE 

01 

DTP • 

•10000000000C 01 

DTG ■ . lOOOOOOOOOOE 

00 

T8L • 

•lOOOOOOOOOOE* 10 

X(13 • -.99270603051E 

00 

X C2) • *107585271 25E 

01 

X (33 ■ 

•57852776536E 

00 

X[4)« ..49AA1502508E' 

•01 

GMM • 

• 12277A71000E*.! 




INTERMEDIATE VALUES 

11 

GOOD STEPS TAKEN # 

0 

BAO STEPS TAKEN 

T ■ .50000000000E 

01 

FPT • .10000000000E 

01 

DTP ■ 

•lOOOOOOOOOOE 

01 

DTG ■ .lOOOOOOOOOOE 

00 

T8L ■ 

•10000000000E*tO 

X (13 • ••18006366363E 

00 

Xt2 ] m *66950196427E 

00 

XC3)« 

-•217861QA633E 

00 

XCM- -.76073847667E 00 

GMM ■ 

•'12277471000E-.1 




INTERMEDIATE VALUES 

91 

0880 STEPS TAKEN , 

o BAD STEPS TAKEN 

T • «60000000000E 

01 

FPT • .lOOOOOOOOOOE 

01 

DTP • 

•lOOOOOOOOOOE 

01 

DTG • .lOOOOOOOOOOE 

00 

T8L ■ 

. lOOOOOOOOOOE* 10 

X(l) • ••216215B3708E 

00 

XC2J- **56708820368E 

00 

XC3]« 

•36135607595E 

00 

XtAJ" ••89A0872M36E 

00 

GMM • 

•12277471000E*01 




INTERMEDIATE VALUES 

u 

G88D STEPS TAKEN , 

0 

BAO 8TEPS TAKEN 

T • •70000000000E 

01 

FPT » .10000000000E 

01 

DTP • 

.lOOOOOOOOOOE 

01 

DTG • .lOOOOOOOOOOE 00 

TOl 9 

. lOOOOOOOOOOE* io 

X(l). •.3731869A31AE. 00 

XC2J • ...106380*79995 01 

X 1 33 * 

••53675A03927E 

00 

XC«f •• IMQ7343734E 00 

GMM 9 

•1227747toqoE*Ol 



TABLE A -2. (Concluded) 


INTERMEDIATE VALUES 9 QBBD STEPS TAKEN , 0 BAD STEPS TAKEN 


T • 

.80000000000E 

01 

FPT • 

• 1 OOQOOOOOOOE 

01 

DTP ■ 

• 1 000 jOOOOOOE 01 

DTQ - , 

•10000000000E 

00 

TOL * 

♦10000000000E*10 

xc»« 

•.93001977239E 

00 

XC2J. 

- • 845521 59396E 

00 

X C33 • 

-•4Q885566983E 0 0 

xt4] ■ 

•52215234765E 

00 

'gmm « 

•1227747I000E-01 






intermediate values 

10 

GOOD STEPS TAKEN # 

0 

BAD STEPS TAKEN 

T • 

• 90000000000E 01 

FPT ■ 

.10000000000E 

01 DTP - 

•10000000000E 

01 

DTQ ■ . 10000000000E 

00 

TOL » •1000000000pE-lO 

xtn- 

-•97913676968E 00 

X[2> 

-•29005836245E 

00 X C3) • 

• 337685853 1 4E 

00 

xc*]f • 4041 441 3046£ 

00 

SMM • • 1227747 IOOqE-OI 






INTERMEDIATE VALUES 

19 

□ 00D STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T - 

•lOOOOOOOOOOE 

02 

FPT • 

.lOOOOOOOOOOE 01 

DTP ■ 

•lOOOOOOOOOOE 

01 

DTG ■ .10000000000E 

00 

TOL • 

•lOOOOOOOOOQE-lO 

xcn- 

-•16592064157E 

00 

XC2l« 

-.3367790P506E 00 

X E 3 3 • 

• 16497137U6E 

01 

X[41« -.20392918830E 

00 

GMM • 

•1 227747 lOOOE-Ol 





INTERMEDIATE VALUES 

31 

GOOD STEPS TAKEN , 

1 

BAD STEPS TAKEN 

T * 

.11000000000E 

02 

FPT • 

.lOOOOOOOOOOE 01 

DTP ■ 

•lOOOOOOOOOOE 

01 

DTG ■ •10000000000E 

00 

TOL » 

•10000000000E*10 

xcn« 

•89885377499E 

00 

X C2I • 

•51450563383E-01 

X C3] ■ 

•681 338641 68E 

00 

XC43- ••12386168192E 

00 

GMM • 

•12277471000E-01 






FINAL 

VALUES 303G00D ! 

STEPS TAKEN 6BAD STEPS TAKEN 


T 

•111243A0337E 

02 

FPT • 

•10000000000E 01 

dtp » 

•10000000000E 

01 

DTG f tlOOOOOOOOOOE 

00 

TBL • 

•lOOOOOOOOOOE-iO 

xcn« 

•9939997551 4E 

00 

X C2] ■ 

• • 897765 352 97E- 06 

X C33 • 

-•146A4710143E- 

-03 

Xt*JT ••20317701066E 01 

GMM ■ 

•12277471000E-01 



TABLE A -3. TEST PROBLEM B12 


^ASSIGN 

S*MT0,SI CR/B0 MT1,L8 UP* 

kREWlND 

MT1 • 


kFGRTRAN 

B8/L8 

• 

1 


dimension xc2i 

2 


DIMENSION ALPHC131,8ETAC13,12J,CHC131 

3 


CO m M0N ALPH, BETA, CH 

4 

O 

fXJ 

Read 11, T I, OTGI , T0L 

5 


READ 1 1, X Cl 1 , X [23 

6 


READ U,TF,FPT1,DTP1 

7 


READ 2025, XT 

8 

2025 

F0OMAT [14] 

9 

C 

CONSTANTS F9R INTEGRATION SUBROUTINE 

10 


08 60 1.1,13 

11 


D0 50 J« 1 , 1 2 

12 

50 

BETA Cl, J1 *0. 

13 


AL°H C 1 1 «0 • 

14 

60 

chcii »o. 

15 


CHC61 -34. 7105. 

16 


CH C 7 3 -9.735. 

17 


CH[g] -CH C71 

18 


CH [91 -9.7280* 

19 


ChCIO] »CH [95 

20 


CHC121 *41 ./8A0. 

21 


CHC133 «CH [ 1 21 

P 2 


Alph C21 -2*727. 

23 


ALPH C31 -1 .79* 

24 


ALPH C43 ■ 1 .76 • 

25 


ALPHC51.5./12. 

26 


alph C63 ■ • 5 

27 


ALPHC73.5./6* 

28 


ALPHC83.1./6* 

29 


ALPH C91 »2 • 73 • 

30 


ALPHC10] «l./3. 

31 


ALPH Cl 1 1 ■ 1 • 

32 


Alph C 1 31 » 1 • 

33 


BETA C2, 11 ■ 2 • /27 • 

34 


BETA C3, 1 1 3 1 ./36. 

35 


BETA C4, 11 -1./24. 

36 


BETA C5, 11 «5./l2. 

37 


BETA C6, 11 * .05 

38 


BETA C7, 11 .-25./108. 

39 


BETA C8, 11 -31 • /300 • 

40 


BETA [9, 11 *2. 

4 1 


BETA CIO, 11 --91./108. 

42 


BETA Cl 1, 11 -2383.74100. 

43 


BETA C12, 11 -3./205. 

44 


BETA Cl 3, 11 •- 1777. 74100. 

45 


BETA C3, 21 -1./12. 

46 


BETA C4, 31 -1./8. 

47 


BETA C5, 31 --25. /16. 

48 


BETAC5,43 «-BETAC5,31 

49 


BETA C6, 41 ■ *25 

50 


BETA C7, 41 -125. /108. 


64 



TABLE A -3. (Continued) 


51 


BETAC9,4l *-53. /6, 

52 


BETAC10,«O»23./l6*. 

53 


BETA Cl 1, 4] *-341. /16*. 

54 


BETAC13,43 *BETAC11,43 

55 


BETAC6,53*«2 

56 


BETA 07,53 *-65. /27. 

57 


BETAC8, 53*61. /225. 

58 


BETA 09, 53 *704. /45. 

59 


BETAtlO,53 *-976. /135. 

60 


BETAC11, 53*4496. /1025. 

61 


BETAC13,53 *BETAC11,53 

62 


BETA 07,63 *125. /5*. 

63 


BETA 08,63 *-2*/9» 

64 


BETA 09,63 *-107. /it. 

65 


BETAC10, 63*311. /B4« 

66 


BETACll,63»-301./82* 

67 


BETA 012, 63 --6./41 . 

68 


BETA 013,63 *-289. /82» 

69 


BETA 08,73 »13./900. 

70 


BETAC9, 73*67. /90. 

71 


BETA 010,73 *-19. /60. 

72 


BETA 011,73 *2133. /4100. 

73 


BETA 012,73 .-3./205. 

74 


BETA 013,73 *2193. /4100. 

75 


BETA 09, 83 *3. 

76 


BETAOIO, 83*17. /6. 

77 


BETAC11, 83*45. /82. 

78 


BETA 012,83 --3./41. 

79 


BETAC13, 83*51. /82. 

80 


BETAOIO, 93*-l. /IE. 

81 


BETACU, 93 *45. /164. 

82 


BETAC12,93*3./41. 

83 


BETA 013,93 »33./164. 

84 


BETAOll, 103*18. /41. 

85 


BETAC12,103 *6./41. 

86 


BETA 013, 103 *12. /41 1 

87 


BET ACl3,123 *1* 

88 


NS»0 

89 


NR * 0 

90 


NST *0 

91 


NRT *0 

92 


PRINT 800 

93 

800 

F0RM AT C1H1,/,51X, 14HINITIAL VALUES3 

94 


Call PRINT CTI,X, FPT1 ,DTP1,DTGI*T0I 

95 


T*T I 

96 


STEP*PPT1 

97 


dtg-fpti-t 

98 


IF CABS CDTG3 • ABS CDTG I) 3 6,6,7 

99 

7 

dtg*dtgi 

100 

6 

NSF»0 

101 


NRF»0 

102 


IFCABSCTF-TI3 -ABS CFPT1-TI3 3 112,121*121 

103 

112 

step-tf 

104 


dtg*dtgi 
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TABLE A -3. (Continued) 

■ 105 121 CALL INTEGRCT,STEP#DTG # T8l#X, 2*KT #NSF*NRF! 

• 136 IFCTF-STEP! 161,151/161 

■ 107 161 T.STEP 

■ 1C8 STEP»T+DTP1 

■ 109 IFfABStOTG! -ABS CDTP1! ! 8*8/9 

■ 110 9 DTG-DTPl 

■ 111 8 IF CA3SCTF-T! -ADSCDTP1!! 132*143*143 

• 112 132 STEP'TF 

■ 1 1 3 IF CABS t DTG! -ABSCTF-Tl! 143/143* 2024 

» 114 2024 DTG.TF-T 

• 115 143 PRJMT 190*NSF*NRF 

• 116 190 F0RMAT tlH *////48X*19H INTERMEDIATE VALUES*3x* I4>1X*19HG88D STEPS T 

• 117 1AKFN / /1X,I4,1X/15HBAD STEPS TAKEN! 

» 118 CALL PRINT CT ,X, FPT1 * DTPi *DTG I , T8L1 

• 119 NS*NS+NSF 

• 1?0 NR.NR+NRF 

• 121 GO TO 121 

• 1?2 151 NS»NS+NSF 

■ 123 NR»NR+NRF 

• 12* NST»NST*NS 

• 125 NRT-NRT+NR 

• 126 PRINT 840/ NST * NRT 

• 127 840 FORMAT C1H0/50X/ 14HFINAL VALUES *I4*16HG88D STEPS TAKEN* J4* 

• 128 115RBAD STEPS TAKEN! 

» 129 CALL PRINTCTF/x* FPTl *DTP1 aDTGI#T8L1 

» 130 GO T6 20 

• 131 H F0RMATr3E18.il! 

» 132 END s •• 


COMMON ALLOCATION 

777*6 ALPH 77256 BETA 77224 CH 

PROGRAM allocation 

00011 X 00015 KT 00016 1 

00020 NS 00021 NR 00022 NST 

00024 NSF 00025 NKF 00026 Tl 

00032 T8L 00034 TF 00C36 FPTl 

00042 T 00044 STEP 00C46 DTG 

subprograms required 
print abs integr 

THE END 


00017 J 
00023 NRT 
00030 DTG I 
00040 DTPI 
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TABLE A -3. (Continued) 


1 

2 

3 

4 

5 

6 

7 

8 
3 


10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
23 


20 


240 

260 

210 

220 

230 


SUBROUTINE I NTEGR CT I # T j DTSj TOL ,X,N#KT#NS,NR) 

DIMENSION XC 2 3iXEC2 3 

nt«o 

NS*0 

NReO 

dtg*dts 

TO-TI 

XEC13 «X [13 
XEC23 *X C23 
STEP-TO+DTG 

Call RK713CT0,STEP # DTG,T0L ,X/ 2 , MS,MR,XE/2 3 

TO*=STEP 

NS«MS+NS 

NRaNR+MR 

dts«dtg 

nt*nt+ms+mr 

IFCSTEP-T3240/230/240 

I F CABS C0TG3 -ABSCT- STEP] 3 210, 210/ 260 

DTG«T-STEP 

IF CNT-KT3 20/220/220 

T “STEP 

RETURN 

END 


PROGRAM ALLOCATION 


DUMMY X 

00030 XE 

00034 INTEGR 

00035 NT 

dummy ns 

DUMMY NR 

00036 MS 

00037 MR 

DUMMY KT 

00040 DTG 

dummy dts 

00042 TO 

DUMMY T I 

00044 STEP 

dummy TOL 

DUMMY T 


SUBPROGRAMS REQUIRED 

RK713 ABS 

the END 
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TABLE A -3. (Continued) 


- 1 SUBROUTINE RK713 CTI,TF,DT , TOL ,X, N, NS*NR,XE,M) 

• 2 c seventh order rungc-kutta integration with stepsize control 

« 3 C TF can be greater than Tl 9R LESS THAN TI AND RK713 WILL WORK 

■ AC NS IS ThC NUMBER OF 3UCCESFULL STEPS TAKEN 

• 5 C NR IS THE NUMBER OF REJECTED STEPS TAKEN 

• 6 C N IS TH E NUMBER OF DIFFERENTIAL EQUATIONS 

• 7 c kt is max \um3fr bf iterations 

• 8 C ARRAY F STORES THE 13 EVALUATIONS PF THE DIFFERENTIAL EQUATIONS 

■ 9 C SUBSCRIPTS FOR ALPHA, BETA, AnD CH ARE *1 GREATER THAN FEHLBERGS 

■ 10 C F[P] IN FEHLBERGS REPORT IS IN FCl,J3 

• 11 C F [ 1 1 IS IN FCl + l.JI 

■ 12 C FEHLBERGS REPORT REFERENCED IS NASA TR R-287 

• 13 c parameters for deq subroutine must be stored in common 

» i4 c dimensions must agree with number of differential equations and 

» is c number of constants in the particular fehlbeRg formula used 

• 16 DIMENSION F C 1 3, 2 ],XDUM[ 2 I, TEC 2 3 ,ALPHtl33, 

• 1 / 1BETA Cl 3, 1 2] , X t 2 I , CH C 1 33 , XE C2 3 

■ 13 COMMON ALPH,BETA,CH 

• 19 T* T I 

• 2 0 NS«0 

» 21 NR ■ 0 

• 22 20 CALL DEQ CX,T,TE3 

• 23 D0 30 I ■ 1 , N 

• 24 30 Fcl, 13 »TECn 

• 25 D0 70 K»2, 1 3 

• 26 DO 40 I ■ 1 , N 

■ 27 40 XDUMCIJ.XCI) 

■ 28 NN«K-1 

• 2S 06 50 1*1, N 

• 30 06 50 J* 1 , NN 

• 31 50 XDUMCn.XDUMCn+OT»BETACK,jD»FCj,I] 

• 32 TDUM*T+ALPH[K3 *DT 

• 33 CALL DEQ CXDUM, TDUM, TE3 

« 34 06 60 I * 1 , N 

» 35 60 F [K, 13 «TE [ 13 

• 36 70 CONTINUE 

■ 37 ER-0. 

. 38 C M IS AN INPUT VALUE WHICH DETERMINES THE NUMBER OF VARIABLES USED IN 

■ 39 C THE ERROR CONTROL LOOP 

• 40 C XE IS AN INPUT VECTOR W,Th DIMENSION M WHICH IS USED TO NORMALIZE 

« 41 C THE TRUNCATION ERROR COMPUTATIONS IN THE ERROR CONTROL LOOP 

• 42 D0 120 I * 1 , M 

■ 43 140 TEm-DT*CFCl,n+FCll,Il-FCl2,n-FC13,in*Al,/sA0./XECn 

« 44 IF CABSCTECI33-ER3 120,120,130 

■ 45 130 ER-A3SCTECI33 

« 46 120 CONTINUE 

» 47 DT1*DT 

« 48 AK» • 8 

• 49 IFCER3141, 142,141 

■ 50 142 DT-10.*DT1 

• 51 GO TO 150 

• 52 141 Dt» rSQRT C SORT C SORT CT8L/ER3 3 3) 

■ 53 DT«AK*DT*0T1 
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TABLE A -3. (Continued) 


54 


IF CER -T0L3 150/150/180 

55 

150 

TF-T+DT1 

56 


08 90 I ■ 1 , N 

o7 


D8 90 L-l/13 

*8 

90 

XCI3 -X CI3+DT1-CHCL3-FCL/I3 

59 


NS-NS+1 

60 


00 T 0 230 

6l 

180 

NR.NR+1 

62 


tf-t 

63 

230 

return 

64 


END 


COWtfo- ALLPC/kTISN 


777*6 

AUPH 

77256 

BETA 

77224 

CH 



PROGRAM 

ALLOCATION! 






0 C C 3 7 

r 

00123 

XDUM 

00127 

TE 

DUMMY 

X 

DU MW V 

XL 

DUMMY 

NS 

Dummy 

NR 

00133 

I 

D U M Y 

N 

00134 

K 

00135 

NN 

00136 

J 

OU Vf/ Y 

M 

00137 

L 

OC 1 4 0 

R<713 

00142 

T 

DUM m Y 

T I 

DUMMY 

DT 

00144 

TDUM 

00146 

ER 

00150 

DTI 

00152 

AK 

dummy 

TOL 

DUMMY 

TF 


subprograms required 

DEG ABS SQRT 

THE FnC 


• 1 Su n RQuT I N'E PRIM [T/X/ FPT/DTP/DTG/T0L3 

« ? DIMENSION X [23 

« ? PRINT l/T/FPT/DTP/DTG/TRL 

■ A 1/XC13/XC23 

• t 1 F8RMAT [1H0/5HT »F. 13*11/ 2X/ 5HPPT ■ E 1 8 . 11 / 2X/ 5HDTP -E18.il/2X/ 

■ 6 1 5hQT S -E18.il/2X/5HTCL -E18.il 

« 7 2///6H XC13-E18.il/2Xi5HXC23-E18.113 

■ 8 RETURN 

■ 9 End 


PROGRAM ALL8CATI8N 


DUM^y x 
OUM V Y OTP 
THE EnD 


OOOIA PRINT 
DUMMY OTG 


dummy t 
dummy T8L 


dummy fpt 
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TABLE A -3. (Continued) 


« 1 SUBROUTINE DEQCX,T,DX] 

« 2 DIMENSION X [23 » DX[23 

■ 3 DXCn«2.*Cl«-XC233*XClJ 

■ 4 DXC2]»-XC23*Cl.-XCt33 

» 5 RETURN 

- 6 END 


PRqGRA^ allocation 

DUMVY X DUMMY DX 00006 DEQ 

THE END 


^ASSIGN BI'MTl. 
aREWIND MT 1 • 
AF0RTL6AD B I U • 
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TABLE A -3 


(Continued) 


NAME ENTRY e»ISl\ last SIZE/10 n M ^0N 


♦PROGRAM 03507 "3477 10024 2262 17063 


Base 


INITIAL VALUES 


T . 

.OOOOOOOO^OOE 

00 

FPT ■ 

.10000000000E 

11 

DTP • 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

01 

DTQ 

• .10000000000E 00 

T8L * 

•10000000000E oo 

xtl)» 

* 1 OOOOOOO^OOE 

01 

x C21 ■ 

.300000000002 

01 












FINAL 

VALUES 4 1 GOOD 

STEPS 

TAKEN 12BAD STEPS 

TAKEN 


T . 

UJ 

o 

o 

a 

o 

o 

o 

o 

o 

o 

o 

ru 

02 

FPT . 

.ICOOOOOOOOOE 

11 

dtp ■ 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

01 

DTG 

• .10000000000E 00 

tol * 

•lOOOOOQOOOQE 00 

xc». 

• 1 1 8324 324 1 6e 

01 

x C 23 - 

.22060344093E 

00 










<1 

to 


TABLE A -3. (Continued) 


INITIAL VALUES 


T . 

xcn ■ 

•COO GO 00 3' OOF 

. icooccoo^doe 

^ 1 

fpt - 
xc?i . 

. loocononoooF 
•30000000000F 

11 

31 

OTP ■ 

•1COOOOOOOOOE 

01 

OTG 

• .lOOOOOOOOOOE 

00 

TOl • 

•IOOOOOOOOOOE-Oi 







FI SAL 

VALUES 32GOOD 

STFPS 

taken isbad steps 

TAKEN 


T • 

xcn» 

.2ococooo' % 30e 

•65544148 °74E 

3? 

.33 

tpt « 
y C2i ■ 

.10003000300E 
•1853463 3 031E 

U 

00 

DTP « 

.lOOOOOOOOOOE 

01 

DTG 

• .lOOOOOOOOOOE 

00 

tol ■ 

•lOOOOOOOOOOE-Ol 








INITIAL V, LUES 




T . 

♦CCOOOCCO^OOE 

00 

FPT • 

.lOOOOOOOOOOE 

11 

DTP - . 1 uOOOOOOOOOE 01 

DTG 

• . 10000000000E 00 TOL • 

•IOOOOOOOOOOE.02 

xcn. 

♦ lOOOOOOOOOOE 

31 

X C?) ■ 

•3000no03000E 

01 










FINAL VALoFS 32G000 

STEPS 

TAKEN UBAD STEPS taken 


T 

• ?3COCCOO A OCE 

0? 

FPT « 

•10000000000E 

11 

otp • .iooooooooooe 01 

DTG 

• .lOOOOOOOOOOE 00 TOl ■ 

•iOOOOOOOOOOE-O? 

xcn. 

.67265 lOS^lE 

00 

x C?' ■ 

• 1P59745R435E 

00 






TAB LE A -3 . ( C ontiriued) 


INITIAL VALUES 


T • 

.OOOOOCOOOOOE 

00 

FPT « 

•10000000000E 

11 

DTP ■ 

.10000000000E 

01 

DTG 

• . 10000000000E 00 

TBL • 

tl0000O00000E-03 

X[1J. 

•10000000000E 

01 

Xt2} • 

.300C0000000E 

01 














FINAL 

VAUUES 38G06D 

STEPS 

Taken iebaD steps 

TAKEN 


T . 

•20000000000E 

02 

FPT ■ 

♦10000000000E 

11 

DTP ■ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

m 

01 

DTG 

■ •10000000000E 00 

T«L ? 

•igoooooooooe-03 

XtU* 

•fe754C936703E 

00 

XI2] ■ 

tl860^35A7if8E 

00 









IMITUU VALUES 


T • 

.OOOOOOOOOOOE 

00 

FPT ■ 

. 10000000000E 

11 

DTP ■ 

.10000000000E 01 

DTG 

• , *10000000000E 00 

T9U • 

• lOOOOOOOOOQEwQA 

xcn« 

aooocooooooE 

01 

X C2} » 

•30000000000E 

01 













final 

VALUES A7G0B0 

STEPS 

taken ubad STEPS ' 

TAKEN 


T 

UJ 

o 

o 

c 

o 

o 

o 

o 

o 

o 

o 

cm 

02 

FpT - 

»10000000000E 

11 

DTP * 

•10000000000E 01 

DTG 

» flOOOOOOOOOOE 00 

T»l ■ 

tlOOOOOOOOOOE-O* 

xcn. 

.67606171007E 

00 

XC21- 

• 1 860794 17 A3E 

00 









TABLE A -3. (Continued) 


INITIAL VALUES 


T • 

•OOCOCOOOOCOE 

00 

ept . 

•lOOOOOOOOOOE 

11 

OTP * 

•1G0000Q0000E 

01 

DTG 

• .10000000000E 00 

TBL • 

•10000000000E.05 

xtn« 

• 1 OOOCOOO-OOOE 

01 

Xt2] • 

•30000000000E 

01 














FINAL 

VALUES 61G88D 

STEPS 

Taken i7Bad steps 

TAKEN 


T > 

.acoocoooncoE 

02 

FPT > 

*10000000000E 

11 

DTP • 

•lOOOOOOOOOOE 

01 

DTG 

■ .lOOOOOOOOOOE 00 

TOL ■ 

•1000000000QE-05 

X(l). 

•67618129°79E 

OC 

Xt2] - 

•186081A6993E 

00 









INITIAL V, LU r $ 


T « 

•OOOOOOOOOCOE 

00 

FPT • 

•lOOOOOOOOOOE 

11 

DTP • 

•lOOOOOOOOOOE 

01 

DTG 

• .lOOOOOOOOOOE oo 

tol • 

•lOOOOOOOOOOE.O* 

XCl3- 

•lOOOOOOOOOOE 

01 

Xt23 ■ 

•30000000000E 

01 














FINAL 

VALUES 75G88D 

steps 

TAKEN 13BAD STEPS 

TAKEN 


T • 

• 20000000 n CCE 

0? 

FPT . 

•10000000000E 

11 

DTP • 

• 1 DOOOOOOOOOE 

01 

DTG 

• .10000000000E 00 

T8L • 

•10000000000E-0* 

xcn* 

•67618638PC3E 

00 

X C?3 • 

• 1860815801 7E 

00 









TABLE A -3. (Continued) 


INITIAL VALUES 


T • 

.ooooooocnccE 

00 

FPT . 

.lOOOOOOOOOOE 

11 

DTP • 

•lOOOOOOOOOOE 

01 

DTG 

• .lOOOOOOOOOOE 00 

T8l p 

•10000000000E.07 

xti)« 

.lOOOOOOOOOOE 

01 

X C2] • 

.30000000000E 

01 














FINAL 

VALUES 96G00D 

STEPS 

TAKEN 11BAD STEPS 

TAKEN 


T . 

tSOOOOOOO^OOE 

02 

FPT * 

•lOOOOOOOOOOE 

11 

DTP ■ 

•lOOOOOOOOOOE 

01 

DTG 

• tlOOOOOOOOOOE 00 

T9U • 

.10000000000E.07 

xcn« 

•67618750253E 

00 

X [23 • 

•18608160742E 

oo 









INITIAL VALUES 


T • 

•OOOOOOOOCCOE 

00 

FPT • 

•lOOQOOOOOOOE 

11 

DTP • 

il^OOOOOOOOOE 01 

DTG 

• • lOOOOOOOOOOE 00 

T0U • 

•lOOOOOOOOOOE.O* 

X tl) • 

•lOOOOOOOOOOE 

01 

XC21 • 

•30000000000E 

01 













FINAL 

VAU’ES 124G00D 

STEPS 

TAKEN 7BAD STEPS ' 

TAKEN 


T . 

•2000C000000E 

02 

FPT ■ 

•10000000000E 

11 

DTP » 

•lOOOOOOOOOOE 01 

DTG 

• tlOOOOOOOOOOE 00 

T0L ■ 

•lOOOOOOOOOOE.O* 

XCl]. 

.676187587Q5E 

00 

XC2)- 

• 1860816091 IE 

00 









<1 

OS 


TABLE A -3. (Continued) 


INITIAL VALUES 


T 

.COOOCOOOCOOE 

00 

FPT * 

,100000000002 

It 

OTP ■ 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

c, 

o 

• 

01 

0TG 

• *100000000002 

00 

•T0l • 

•10000000000E-09 

xcn« 

•looocooorcoE 

01 

XC2) * 

.30000000000E 

01 















FINAL 

VALUES 162GG0D 

STEPS 

TAKEN 63A0 STEPS 

TAKEN 


T 

• 20000000 n CCE 

02 

FPT • 

•10000000000E 

11 

DTP ■ 

.100000000002 

01 

DTG 

LJ 

O 

O 

o 

o 

o 

o 

o 

o 

o 

o 

• 

00 

t»l • 

•10000000000E-09 

*cn- 

.6761376Q770E 

00 

x C2) • 

.18608160965E 

00 










initial VALUES 


T 

*00000000^002 

00 

FPT • 

. ioooooocoooe 

11 

DTP • 

.100000000002 01 DTG 

* *100000000002 00 

T0l • 

• 10000000000CMO 

X Cl3 - 

.1 0000000^002 


X C 23 * 

,300000000002 

01 

final 

VALUES 217GOOD STEPS 

TAKEN 15BAD STEPS 

TAKEN 


T 

,2oooooooooce 

02 

FPT • 

,100000000002 

11 

DTP • 

.100000000002 01 DTG 

• .10000000000E 00 

T0l • 

.lOOOOOOOOOOE.tO 

X Ct 3 • 

.f,76!8761?79E 

00 

Xt?l- 

.186081609612 

00 







TABLE A -3 


(Continued) 


INITIAL VALUES 


■ ■ 

*- X 

•OOOOCCOO^OCE 00 
. 100000000.CCE 01 

FPT • 
XC23 ■ 

•10000000000E 

•30000000000E 

01 DTP • ,100000000002 

01 

01 

DTG 

■ .IOOOOOOOOOOE 

00 

T8l • *100000000002-10 





INTERMEDIATE values 


GQ90 

STEPS TAKEN $ 

2 

BAD STEPS TAKEN 

T 

XC13- 

.IOOOOOOOOOOE 01 
• 773 A40 1 61 5AE-01 

FPT ■ 
X C2] ■ 

.lOOOOOOCOOOE 
• H644A81571E 

01 DTP « ,100000000002 

01 

01 

DTG 

- .IOOOOOOOOOOE 

00 

TOl • *100000000002-10 





INTERMEDIATE VALUES 

10 

GOOD 

STEPS TAKEN $ 

0 BAD STEPS TAKEN 

T 

xcn« 

•200O0COO0COE 01 
.84977753179E-01 

FPT • 
XC23 ■ 

.10000000000E 

•57795270690E 

01 DTP • .IOOOOOOOOOOE 
00 

01 

DTG 

* .IOOOOOOOOOOE 

00 

TBU ■ ilOOOOOOOOOOE-lO 





INTERMEDIATE VALUES 

9 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

xcl3- 

• 30000000000E 01 
•290891351 73E 00 

FPT . 
X C23 « 

•iooooooooooe 

• 2*92531 7274E 

01 DTP « .IOOOOOOOOOOE 
00 

01 

DTG 

• .IOOOOOOOOOOE 

00 

TBl • *10000000000E-10 





INTERMEDIATE VALUES 

9 

GOOD 

STEPS TAKEN » 

1 

BAD STEPS TAKEN 

T 

X C 1] ■ 

•40000000000E 01 
• 1 A A66020925E 01 

FPT • 
XC2] ■ 

•iooooooooooe 

•18721896503E 

01 DTP • .lC'OOOOOOOOOE 
00 

01 

DTG 

• .IOOOOOOOOOOE 

00 

TBt ■ *100000000002-10 





INTERMEDIATE VALUES 

u 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

X [1 3 ■ 

•50000C00OC0E 01 

• A05HA706*8E 01 

FPT • 
X C2] » 

.iooooooooooe 

• U39A903986E 

01 DTP ■ .IOOOOOOOOOOE 
01 

01 

DTG 

» .IOOOOOOOOOOE 

00 

TBl ■ *10000000000E-lO 





INTERMEDIATE VALUES 

19 

GOOD 

STEPS TAKEN , 

0 

BAD STEPS TAKEN 

T 

X£13- 

•600000000COE 01 
•17561A72735E 00 

FPT - 
X [23 • 

•iooooooooooe 

•2258589A715E 

01 DTP ■ .IOOOOOOOOOOE 
01 

01 

DTG 

■ .IOOOOOOOOOOE 

00 

T8l • *10000000000E-10 





INTERMEDIATE values 

12 

GOOD 

STEPS TAKEN # 

1 

BAD STEPS TAKEN 

T • 

X[13« 

• 70000000000E 01 

•65310426683E-01 

FPT - 
XC23- 

• iooooooooooe 
•90879526320E 

01 DTP - .IOOOOOOOOOOE 
00 

01 

DTG 

■ .IOOOOOOOOOOE 

00 

TOL " .10000000000E-10 


<1 

<1 


<1 

00 


TABLE A -3. (Continued) 



INTERMEDIATE VALUES 8 GOOD STEPS TAKEN < 0 BAD STEPS TAKEN 

T • .800OCC0OOO0E 01 F PT . , 1 OOOOOOOOQOE 01 DTP • .lOOOOOOOOOOE 01 DTG T • 10000000000E 00 TOl • • lOOOOOOOOOOE- iO 

XC13- .1*72?6B2^0SE 00 XC2 ]■ ,36671‘S835?7E 00 


INTERMEDIATE VALUES 8 GOOD STEPS TAKEN , 0 BAD STEPS TAKEN 

T . .90DOCCOO"OOE 01 FPT « . 1 OOCCOOOOOOE 01 DTP • •lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 TOL • . 10000000000E*10 

Xtn- .65059S55881E 00 X [2l - . 1 8757387506E 00 


INTERMEDIATE VALUES 11 GOOD STEPS TAKEN , 0 BAD STEPS TAKEN 

T . .lOOOOOOOOOOE 02 FPT - .lOOOOOOOOOOE 01 DTP ■ .lOOOOOOOOOOE 01 DTG ■ .lOOOOOOOOOOE 00 TOl ■ . lOOOOOOOOOOE.lO 

XCD« • 31 443368°! OE 01 XC2]- .34881916525E 00 


INTERMEDIATE VALUES 20 GOOD STEPS TAKEN , 0 BAD STEPS TAKEN 

T - .HOOCCCOPOOE 02 FPT - .lOOOOOOOOOOE 01 DTP ■ .lOOOOOOOOOOE 01 DTG ■ •lOOOOOOOOOOE 00 T8 l • • 10000000000E*10 

XC13- .90951832589E 00 Xt2l» •29967302638E 01 


INTERMEDIATE VALUES 14 G68D STEPS TAKEN , 0 BAD STEPS TAKEN 

T . • 1 200CC00000E 02 FPT • , lOOOOOOOOOOE 01 DTP ■ •lOOOOOOOOOOE 01 DTG • •lOOOOOOOOOOE 00 TGl ■ • 10000000000E*10 

XC11- .757 1545503 6E -01 X C2l • . 14327141371E 01 


INTERMEDIATE VALUES 9 GOOD STEPS TAKEN , 0 BAD STEPS TAKEN 

T . .13000000°COE 02 r PT . , 1 OOOOOOOOQOE 01 DTP • . 1 OCOOOOOOOOE 01 DTG • • lOOOOOOOOOOE 00 TOl • . 10000000QOOE*lO 

X(l] . .86722147974E-01 Xt2l • . 56555380 1 1 8E 00 


INTERMEDIATE VALUES 7 G50D STEPS TAKEN , 0 BAD STEPS TAKEN 

T . .lACOOCOOOOOE 02 FPT ■ .lOOOOOOOOOOE 01 DTP ■ .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 T8t • . 1000000000QE-10 

X clj - • 301 46933561E 00 X C2l * *245l 25791 98E 00 

INTERMEDIATE VALUES 9 GOOD STEPS TAKEN , 0 BAD STEPS TAKEN 

T - . 15000000000E 02 FPT > .lOOOOOOOOOOE 01 DTP - .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 T&L • . IOOOOOOOOOQEpIO 



TABLE A -3. (Concluded) 


xcn- 

.15034C3a838E 

01 

XC8] ■ 

,189339814$BE 

o 

o 



i 4 






INTERMEDIATE VALUES 

16 

GOOD 

STEPS TAKEN # 0 BAD STEPS TAKEN 

T . 

•16000000000E 

0? 

fpt ■ 

t 10000000000E 

01 DTP ■ tlOOOOOOOOOOE 

01 

DTG 

■ tlOOOOOOOOOOE 00 f T8L « .10000000000E*10 

xcn« 

•395790397^5E 

01 

X C2] • 

• 1 5458979640^ 

01 










INTERMEDIATE values 

80 

GOOD STEPS TAKEN , 

0 bad STEPS taken 

T • 

•17000000000E 

02 

FPT ■ 

.10000000000E 

01 DTP « tlOOOOOOOOOOE 

01 

DTG » .10000000000E 

00 TOL • .10000000000E-iO 

X ClJ « 

.16560785758E 

00 

XI2T 

• 22U580A662E 

01 








intermediate values 

13 

GOOD STEPS TAKEN » 

1 BAD STEPS TAKEN 

T ■ 

X ClJ • 

•laoooooooooE 02 
•6562A78S633E-01 

FPT • 
XC23 ■ 

.10000000000E 

•88886887967E 

01 DTP • • 10O00000000E 

00 

01 

DTG ■ tlOOOOOOOOOOE 

00 TOL ■ *10000000000E»lO 







intermediate values 

8 

GOOD 

STEPS TAKEN , 

0 BAD STEPS TAKEN 

T 

• 1 9000000000E 

02 

FPT • 

.10000000000E 

01 

DTP • 

tlOOOOOOOOOOE 

01 

DTG 

« tlOOOOOOOOOOE 

00 TOL ■ . lOOOOOOOOOOE-10 

xtn- 

.1517AA137A7E 

00 

X t23 • 

.35939068967E 

00 












final 

VALUES 839G00D 

STEPS 

TAKEN 5BAD steps taken 

T . 

•20000000000E 

08 

FPT . 

,10000000000E 

01 

DTP - 

tlOOOOOOOOOOE 

01 

DTG 

■ tlOOOOOOOOOOE 

00 TOL • .10000000000E.10 

xcn - 

• 6761 8761 ?88E 

00 

x cea ■ 

•18608160960E 

00 







TABLE A -4. TEST PROBLEM E22 


**SS!.',N S*HTO,SI CR, Bfi mt1,e8 |_P« 

AEEWi\D MT 1 • 

AF9RTHAN 38, L8. 

■ 1 Ol'-ENSieS X [23 

« 2 DIMENSION ALPHC13}.«3ETAC13»123/CHCJ33 

• 3 COMMON, A(_PH,BEta,Ch 

• 4 20 REAO 11, TI,dTGI,T8L 

« 5 PE AO ll,XC13,Xr23 

• 6 READ ll,TF,FPTl,OT»’l 

• 7 RE *0 2025, <T 

• 8 2025 F 0 DV AT f 1 41 

• 9 C C8'3TA\TS FOR INTE3RAT19N sUBRSUTINE 

• 10 08 60 1*1,13 

• 11 Of) 5C J«l,12 

« 12 50 BETA Cl, Jl *C. 

• 13 Al-HCHO. 

■ 14 60 CM [ 1 3 * 0 • 

• 1 5 CH r 63 *74 ./105. 

• 16 Ch f 73 *9 . /35 • 

• 1 7 Cm [83 *CnC73 

« 18 CM r 93 * 9 • /2?0 • 

• 19 ChC 103*CMC o 3 

• 2C CHC123 *41 ./8A0. 

• 21 CH C 133 *CH C 1 23 

• 22 Al°HC23»2./27. 

« ~3 Al^H C33 *1 ./9» 

» 24 ALPH [43 = 1 ./6* 

• ?b A L ph [53 *5./12* 

• 26 Al°H[63».5 

» 27 AL°m [73 *5./6. 

« 28 A l ph[S3 =1 ./6« 

• 29 AlPm [93 *2./3* 

• 30 Ai_ph [1C3 * 1 «/3* 

» 31 AL D H [113*1. 

• 32 Al d h [133*1* 

« 33 BE T A C2, 13 *2./27. 

■ 34 BETAC3, 13 *1 ./36. 

« 95 BETA [4, 13 *1 ./24. 

• 36 BETA [5, 13 »5./l?. 

• 37 BETA [6, 13 * .05 

• 98 BE T A C7, 13 «-25./10-? . 

« 39 3ETA [g, 13 *31 ./900. 

. 40 BETA(9,13«2. 

• 4 1 BETA [1C, 13 *"91 • / 1 0 8 • 

« 42 BETA [1 1, 13 *2383./ -100. 

« 43 BETA [12, 13 *3*/205- 

• 44 BETA [13, 13 *-1777. /4100. 

• 45 BETAC3,23 *1./12. 

• 46 BETA[4,33*l./8. 

« 47 BETA C5,33 --25./16. 

« 48 BETAC5,43*-BETa[5.33 

• 49 Beta £6,43 • .25 

• 50 BETAC7,43 -125./108. 
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TABLE A -4. (Continued) 


51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 
67 
6 8 

69 

70 

71 
73 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 
8 3 
8% 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 


BETA [9/41 --53. /6. ' 

BETAC10, 43*23. /108. 

BETACU/41 --341 ./164. 
3ETAC13,43*BETA[ll/43 
BETAC6/51 * • 2 
BETA [7/53 «-65./27. 

BETA [8/53 *61 ./225. 

BETA [9/53 *704. /45. 

BETA [10/53 --976. /135. 

BETAC11, 53*4496. /1025. 

BETA [13,53 »SETA(11,53 
BETA [7/63 *125. /54. 

BETA [8/63 * ”2 • /9 » 

BETA [9/ 63 *-107. /9. 

BETA [10/63 *311 ./54 . 

BETA [1 1, 63 *-301 ,/82. 

BETA [12, 63 *-6»/41. 

BETA [13/63 *-289. /82» 

BETA [8/ 73 *13./900« 

BETA (9/73 *67 • /90 • 

BETA [10, 73 *-19. /6C. 

BETAdl, 73 *2133. /4100. 
BETA(12,73*-3./205. 

BETA [13,73 *2193. /4i00* 

BETA [9/81 *3. 

BETA[10,8] *17. /6. 

BETA (1 1,83 *45. /82. 

BETA [12,81 *-3./41 . 

BETA (13, 81 *51 */82. 

BETA (10,93 *-l ./12. 

BETA [1 1/93 -45./164 . 

BE t A (12/ 93 * 3 • /4 1 • 

BETA (13/93 *33. /164 . 

BETA (11, 101 *18. /41 . 
beta ( 12 , 101 -6./41 . 

BETA [13, 103 *12. /41 . 

BETA (13, 123 *1. 

NS = 0 
NR*0 
NST*0 
NRT » 0 
PRINT 800 

800 F6PNAT (1H1///51X, 14HINIT1AL VALUES3 

CALL PRINT [TI/X, FPT1 / DTP! , DTG I , T813 
T.TI 

STFP*FPT1 

DTB*FPT1-T 

IF [ABS (DTG3 -ABS (DTG 13 3 6/ 6/ 7 
7 DTG* DTG I 

6 NSF *0 

NRF*0 

IF CABS (TF-TI3 - ABS [8-PT1-T 13 3 112/12li 121 

112 STFP*TF 
DTG*DTGI 



TAB LE A -4 . ( C ontinued) 


• 

105 

121 

CALL INTEGRCT/STEP/DTG/T8t./X, g,«T /NSF/NKF) 

• 

106 


IFCTF-STEP)161/151/161 ... 


107 

161 

T.STEP 


108 


step«t+dtpi 

m 

109 


IFCABSCDTG)-ABSCDTPl))8/8/9 


110 

9 

DTG-DTP1 


111 

8 

IF tABS CTF-T) -ABS CDTP1) ) 132/143*143 


112 

132 

STEP-TF 


113 


IF CABS CDTGJ -ABS CTF-T) ) 143/1*3/ 2024 


114 

2024 

DTG-TF-T 


115 

143 

PRINT 1 90/ NSF / NRF 


116 

190 

FORMAT C1H //// > 48X / 19H I NTERMEDI ATE VALUES/ 3X/ 14# 1 X/ 19HG090 STEPS T 


117 


1AKFN / /1X/I4/1X/15HBAD STEPS TAKEN! 


118 


call print it , x , fpti/Dtpi/Dtgi/Tod 


119 


NS-NS+NSF 


120 


NR-NR+NRF 


121 


GO TO 121 


122 

151 

NS.NS+NSF 


123 


NR-NR+NRF 


124 


NST.NST+NS 


125 


NRT-NRT+NR 


126 


PRINT 840/ NST/ NRT 


127 

840 

FORMAT C1H0/50X/ 14HFINAL VALUES /I4/16HG00D STEPS TAKEN/ I*/ 


128 


115HBAD STEPS TAKEN) 


129 


CALL PRINTCTF/X/ FPTl / DTPl / DTGl / TOL) 


130 


GO TO 20 


131 

11 

F0RMATC3El8.il) 


132 


END 


C0MM0N ALLOCAT ION 

77746 ALPH 77256 BETA 77224 CH 

PROGRAM allocation 

OOOll X 00015 KT 00016 I 00017 J 

00020 NS 00021 NR 00022 NST 00023 NRT 

0002* NSF 00025 NRF 00026 TI 00030 OTQI 

00032 TOC 00034 TP 00036 FPT1 000*0 DTP1 

00042 T 00044 STEP 00046 DTG 

SUBPROGRAMS REQUIRED 

print abs INTEGR 

the end 
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TABLE A -4. (Continued) 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


20 

21 

22 

23 


20 


240 

260 

210 

220 

230 


SUBROUTINE I NTEGR CT I t T, DTS* TOL *X*N,KT*NS*NR) 

DIMENSION Xt 2 3*XEC2 3 

NT = 0 
NS = 0 
NR* 0 

dtg-dts 
TO = TI 
XE C13 »1. 

XEC23 »1. 

STEP=TO+DTG 

CALL RK713CT0/STEP # DTG,T8L ,X* 2 t MS*MR,XE*2 3 

TO-STEP 

NS-MS+NS 

NR *NR + MR 

DTS»DTG 

nt-nt+ms+mr 

IF CSTEP-T3 240*230*240 

IFCABSCDTG3 -ABSCT- S fEP3 3 210*210*260 

DTG * T -STEP 

IFCNT-KT3 20*220*220 

T» STEP 

RETURN 

END 


PROGRAM ALL0CAT I 8N 


dummy x 

00030 XE 

00034 INTEGR 

00035 NT 

dummy ns 

DUMMY NR 

C0036 MS 

00037 MR 

dummy kt 

00040 DTG 

dummy dts 

00042 TO 

DUMMY T 1 

00044 STEP 

dummy tbl 

dummy T 


subprograms reguired 
RK713 ABS 

the end 
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TABLE A -4. (Continued) 


• 1 SUBROUTINE RK713 IT I / TF/ DT ,T0L ,X,N, NS/NR/XE/M) 

■ 2 c seventh order runge-kutta integration with stepsize control 

■ 3 C TF CAN BE GREATER THAN Tl 9R LESS THAN TI AND RK713 WILL WORK 

■ 4 C NS IS THE NUMBER OF SUCCESFULL STEPS TAKEN 

« 5 C NR IS THE NUMBER 0F REJECTED STEPS TAKEN 

• 6 C N IS THE NUMBER BF DIFFERENTIAL EQUATIONS 

■ 7 C KT IS MAX NUMBFR OF ITERATIONS 

■ 8 C ARRAY F STORES THE 13 EVALUATIONS OF THE DIFFERENTIAL EQUATIONS 

« 9 C SUBSCRIPTS FOR ALPHA, BETA/ AND CH ARE +1 GREATER THAN FEHLBERGS 


10 

C 


F C01 in FEHLBERGS REPORT IS IN FCl,J3 


11 

c 


Ren is in Fn + i/ji 


12 

c 


FEHLBERGS REPORT REFERENCED IS NASA TR R-287 

13 

c 


parameters for deq subroutine must be 

STORED IN COMMON 

H 

c 


DIMENSIONS must AGREE WITH number OF DIFFERENTIAL EQUATIONS and 

15 

c 


number of constants in the particular 

FEHLBERG FORMULA USED 

16 



dimension F [13/ 2 1/XDUMC 2 3/TEt 2 1 

/ALPHtl33/ 

17 



1BETA C13/12I ,XC 2 3 ,CHtl33 ,XEC2 3 


18 



COMMON ALPH/BETA/CH 


19 



T.TI 


20 



NS*0 


21 



nr»o 


22 


20 

CALL DEQ CX/T/TE3 


23 



DO 30 I * 1 , N 


24 


30 

F [ 1 / 1 3 * TE C 13 


25 



DO 70 K-2/13 


26 



DO 40 I « 1 / N 


27 


40 

XDUM C 1 3 -X CI3 


28 



NN«K-1 


29 



DO 50 I.l/N 


30 



DO 50 Js 1 / NN 


31 


50 

XDUMCI3 «XDUMCI3+DT*BETACK/J3*FCJ, I) 


32 



TDUM-T + ALPHCK3 *DT 


33 



CALL DEQ CXDUM/TDUM/TE3 


34 



DO 60 I ■ 1, N 


35 


60 

FCK/I3-TECI3 


36 


70 

CONTINUE 


37 



ER*0» 


38 

c 

M 

IS AN INPUT VALUE WHICH dete Rm ines the 

number of variables used II 

39 

c 

THE ERROR CONTROL LOOP 


40 

c 

XE 

IS AN INPUT VECTOR WITH DIMENSION M WHJCH IS USED TO NORMALIZE 

41 

c 

THE TRUNCATION ERROR COMPUTATIONS in the 

error CONTROL loop 

42 



DO 120 I * 1 / M 


43 

140 

TEC II «DT* CFClz IJ+Ftll>I3-Fll2# I3«Ftl3, 

133*41, /840./XECI3 

44 



IF (ABSCTECI33-ER3 120/120/130 


45 


130 

ERo ABS CTE 1 1 3 3 


46 


120 

CONTINUE 


47 



dti«dt 


48 



AK»»8 


49 



IFCER3 141/142/141 


50 

142 

DT»10.#DT1 


51 



GO TO 150 


52 


141 

Dt» tSQRTtSQRTtSQRTCTOL/'ERIHl 


53 



dt«ak*dt«dti 
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TABLE A -4. (Continued) 


■5* IF CER -T0L1 ISOaITOaISO 

55 150 TF-T+DT1 

56 Da 90 I-1 aN 

57 06 90 L* 1a 13 

58 90 XCn»XCn+DTl*CHCL3*rCLiI3 


59 


NS«NS+1 

60 


GO TO 230 

61 

180 

NRmNR+1 

62 


TF-T 

63 

230 

RETURN 

64 


END 


COMMON allocation 


77746 

ALPH 

77256 

BETA 

program 

ALLOCATION 


OOC 37 

F 

00123 

XDUM 

dum^y 

XE 

DUMMY 

NS 

DUM^Y 

N 

00134 

K 

Dummy 

M 

00137 

L 

DUm^Y 

T I 

DUMMY 

DT 

00150 

DTI 

00152 

AK 


77224 CH 


DC 1 27 

TE 

DUMMY 

X 

DUMMY 

NR 

00133 

I 

00135 

NN 

00136 

J 

0C140 

RK713 

00142 

T 

Ot'l 44 

TDUM 

00146 

ER 

DUMMY 

TOL 

DUMMY 

TF 


SUBPR°ORAMS required 

DEQ ABS SORT 

THE enD 


« 1 

■ 2 

• 3 

> 4 

» 5 

« 6 

» 7 

« 8 

■ 9 


SUBROUTINE PRINT CTaXa FPTaDTPaDTGaTOL) 

DIMENSION X C23 
PRINT IaTaFPTaDTPaDTCaTOL 
1aXC13aXC2J 

FORMAT Cl HO/ 5HT « El 8 • 1 1 • 2Xa 5HFpT > El 8 . 1 1 a 2Xa 5HDTP ■E18*lf»2X/ 

15HDT3 ■E18.11a2X,5HT0L »E18*U 
2/ / » 6H XC11 *E18# 11,2X,5HX C23 *E18.111 

return 

end 


program ALLOCATION 


dummy X 
DUMMY DTP 

The end 


00014 PRINT 
DUMMY DT6 


DUMMY T 

dummy tol 


dummy fpt 
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TABLE A -4. (Continued) 


• 1 SUBROUTINE DEQCX/TVDX3 

• 2 DIMENSION XC23, DX C23 

• 3 DXC13-XC23 

• 4 DXC23«C1,-XC13**23#XC23-XC13 

• 5 RETURN 

■ 6 END 


PROGRAM allocation 

DUMMY X DUMMY DX 00006 DEO 

the end 


aasSIGN BI-MJI. 
AREWIND mti • 
aFQRTLDAD biu. 
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TABLE A -4. (Continued) 


NAME entry origin cast size/io common base 

•PROGRAM 03507 03477 10637 2657 17063 


INITIAL VALUES 


T • 

•OOOOOOOOOOOE 

00 

FPT . 

•IOOOOOOOOOOE 

11 

DTP » 

•IOOOOOOOOOOE 01 

DTG 

r • iooooooooooe oo 

TOl • 

•iooooooooooe oo 

xtn. 

•20000000000E 

01 

X C2] • 

•OOOOOOOOOOOE 

00 












FINAL 

VALUES 35G8SD 

STEPS 

TAKEN 13BA0 STEPS TAKEN • 


T . 

•20000000000E 

02 

FPT • 

•10000000000E 

11 

DTP ■ 

.IOOOOOOOOOOE 01 

DTG 

> •IOOOOOOOOOOE 00 

TOL • 

• 100000000** 00 

Kill- 

•19735582007E 

01 

XC23- 

• M827137254E 

00 







00 

<1 



00 

00 


TABLE A -4. (Continued) 


INITIAL VALUES 


T • 

• 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

n 

00 

tpt ■ 

•10000000000E 11 

DTP ■ • 1 OOOOOQOOOOE 01 

,DTG 

■ •10000000000E 00 T&L * 

.lOOOOOOOOOOE«Ot 

1 tl) * 

•20GOOOOOOOOE 

01 

XC2J# 

•OOOOOOOOOOOE 00 










FINAL VALUES 36Q88D ! 

STEPS 

TAKEN 15BAD STEPS TAKEN 


T ■ 

.20000000000E 

02 

FPT ■ 

t 10000000000E 11 

DTP • .10000000000E 01 

DTG 

T tl.OOOOOOOOOOE 00 T6L ■ 

•tOOOOOOOOOOE«Ot 


•20069643660E 

01 

Xt23» 

• *718450490 17E* 01 





; 





INITIAL VALUCS 




T 

•OOOOOOOOOOOE 

00 

FPT • 

flOOOOOOOOOOE 11 

DTP ■ .10000000000E 01 

DTg 

*10000000000E 00 tol • 

. •tooooooooooc*ot 

XtUm 

•20000000000E 

01 

X C2l ■ 

•OOOOOOOOOOOE 00 










FINAL VALUES 37Q08D 

STEPS 

TAKEN 179AD STEPS TAKEN 


T - 

.20000000COOE 

02 

FPT • 

•10000000000E 11 

DTP • .10000000000E 01 

DTQ 

• , 10000000000E 00 T8L ■ 

.1000000000#f.*t 

xcn« 

•20081874022E 

01 

XC23- 

-•41737864121E-01 







TABLE A -4 


(Continued) 


INITIAL values 


T « 

•OOOOOOOOOOOE 

00 

FPT - 

•IOOOOOOOOOOE 

11 

DTP ■ 

•10000000000E 

01 

DTG 

■ .10000000000E 00 T&L * 

•10000000000E-03 

xcn- 

• EOOOOOOOOOOE 

01 

X C23 • 

•OOOOOOOOOOOE 

00 













FINAL 

values a6G99D 

STEPS 

Taken 16BAD STEPS TAKEN 


7 • 

•200QOOOOOOOE 

02 

FPT • 

. 10000000000E 

11 

DTP * 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

01 

DTG 

* * 10OOOOOOOOOE 00 Tfl^ ■ 

•10000000000E-03 

xcn» 

• 2008 1 6A4290E 

01 

X C?] ■ 

-•42197190845E- 

■01 







INITIAL VALUES 


T • 

•OOOOOOOOOOOE 

00 

FPT ■ 

•10000000000E 

11 

DTP ■ 

•iooooooooooe 

01 

DTG 

■ •IOOOOOOOOOOE 00 

tol ■ 

•10000000000E-0^ 

xcn* 

• EOOOOOOOOOOE 

01 

x C2] • 

•OOOOOOOOOOOE 

00 














FINAL 

VALUES 56GO0D 

STEPS 

taken 180AO STEPS TAKEN 


T . 

•20000COOPOOE 

02 

FPT ■ 

•10000000000E 

11 

DTP ■ 

•IOOOOOOOOOOE 

01 

DTG 

* •IOOOOOOOOOOE 00 

TBt ■ 

•iooooooooooe.o* 

xcn« 

•20081537041E 

01 

XC21 • 

..A250A030549E. 

-01 









CO 

o 


TABLE A -4. (Continued) 

INITIAL VAtUES 

T • .OOOOOOOOOOOE 00 FPT ■ « IOOOOOOOOOOE 11 OTP ■ • IOOOOOOOOOOE 01 DT3 • • 10000000000E 00 Tdl • • IOOOOOOOOOOE* 06 

XC13 * * •2QOOOOQO0OOE Oi X C21 « * ' *OOOOOGOOOOOE 00 * 

FINAL VALUES 70GMD STEPS TAKEN 18BAD STEPS TAKEN 

T • *20000000000£ 02 TPT • .IOOOOOOOOOOE 11 OTP » .IOOOOOOOOOOE 01 OTQ • * IOOOOOOOOOOE 00 T8t • • IOOOOOOOOOOE *05 

X Cll • .20081500503E 01 XC2l* • .42508225085E-01 


INITIAL VALUES 

T * .OOOOQOOOOOOE 00 FPT • .IOOOOOOOOOOE 11 OTP • .IOOOOOOOOOOE 01 ,DTQ ?..• IOOOOOOOOOOE 00 T®L * • tOOOOOOOOOQE^OS 

KtU* • 20000000000E 01 Xt2l« .OOOOQOOOOOOE 00 

FINAL VALUES 89G800 STEPS TAKEN 21BAD STEPS TAKEN 

T • .20000000000E 02 FPT ■ .IOOOOOOOOOOE U OTP ■ .IOOOOOOOOOOE 01 OTQ f * IOOOOOOOOOOE 00 T&L • *10000000000E«OS 

X (11 ■ *20081 49788AE Ol KC2)« - 1 A25Q8639103E-01 


TABLE A -4 


(Continued) 


INITIAL VALUES 


T . 

.OOOOOOOOOOOE 

00 

c-PT • 

.IOOOOOQOOOOE 

11 

DTP ■ 

• iOOOOOQOOOOE 01 

DTQ 

• .10000000000E 00 

tol • 

•10000000000E-07 

Xtl3» 

•20000000C00E 

01 

XC23« 

•OOOOOOOOOOOE 

00 













FINAL 

VALUES 112G88D 

STEPS 

TAKEN 18BAD STEPS TAKEN 


T . 

•20000000000E 

02 

FPT * 

•10000000000E 

11 

DTP • 

•10000000000E 01 

DTQ 

* .IOOOOOQOOOOE 00 

tol ■ 

•10000000000E*07 

XCl)* 

•20081*976502 

01 

X C2) ■ 

-.*2508870*142- 

'01 












INITIAL VALUES 




t • 

•OOOOOOOOOOOE 00 

FPT • 

.IOOOOOQOOOOE 11 

DTP • .IOOOOOQOOOOE Oi 

DTG 

• .IOOOOOQOOOOE 00 TOL • 

• lOOOOOOOQOOE.M 

xtn. 

•20000000000E 01 

X C2] » 

.OOOOOOOOOOOE 00 

FINAL VALUES 1*%G88D 

STEPS 

TAKEN 16BAD STEPS TAKEN 


T « 

•20000000000E 02 

FPT . 

•10000000000E 11 

OTP • .10000000000E 01 

DTQ 

• .IOOOOOQOOOOE 00 T0L • 

• IOOOOOOOOOoE.M 

X(l]« 

.20081*97620E 01 

X C2) « 

-.42508880677E-01 






zo 

to 


TABLE A -4. (Continued) 


INITIAL VALUES 


T . 

•OOOOOOOOOOOE 

00 

FPT • 

. 10000000000 E 

11 

OTP • 

. 10000090000 E 01 

OTfl 

f . 10000000000 E 00 

T 8 U ■ 

.looooooooooe.o* 

utn» 

* •20000000000E 

01 

V C25 ■' 

•OOOOOOOOOOOE 

00 












FINAL 

VALUES 186G90D 

STEPS 

TAKEN 12BAD STEPS 

taken 


T . 

•aooooooooooE 

02 

FPT . 

tlOOQOOOOOOOE 

11 

DTP « 

•10000000000E 01 

OTO 

• • lOCJOOOOOOOE 00 

T9L t 

•ioooooooooqe«os 

xtn* 

.20081497617 E 

01 

X C2) • 

-•4250S884626E' 

-01 








INITIAL VALUES 


T . 

•OOOOOOOOOOOE 

00 

FPT ■ 

.10000000000E 

11 

DTP • 

•10000000000E 01 

OTQ 

■ f , • 10000000000 E 00 

tol • 

.IOOOOOOOOOOE.IO 

Xtn- 

•20000000000E 

01 

XC23« 

•OOOOOOOOOOOE 

00 












final 

VALUES 242G06D 

STEPS 

TAKEN 4BAD STEPS 

taken 


T • 

•20000000000E 

02 

FPT 1 

•10000000000E 

11 

DTP ■ 

• 10000000000 E 01 

OTQ 

• tlOOOOOOOOOOE 00 

T8l • 

.IOOOOOOOQOOE.1. 

xcn- 

•20081497615E 

01 

XC23 ■ 

.t4250R8865l7E- 

-01 









TABLE A -4. (Continued) 


INITIAL VALUES 

T • .OOOOQOOOOOOE 00 FPT • .lOOOOOOOOOOE 01 DTP • . 1 OOOOOOOOOOE 01 DTG ■ .lOOOOOOOOOOE 00 T0L ■ •10000000000E*iO 
X C 13 « .20000000000E 01 Xt2]« .OOOOOOOOOOOE 00 


INTERMEDIATE VALUES 13 GOOD STEPS TAKEN , 1 BAD STEPS TAKEN 

T • .10000000000E 01 FPT ■ • lOOOOOOOOOOE 01 DTP ■ . 1 OOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 TOL • .iOOOOOOOOOOE.tO 

XCU* .15081A42366E 01 X C2) ■ - .78021807484E 00 


INTERMEDIATE VALUES 13 GOOD STEPS TAKEN » 0 BAD STEPS TAKEN 

T • .20000000000E 01 FPT • .lOOOOOOOOOOE 01 DTP • • 1 OOOOOOOOOOE 01 DTG ■ ,100000000002 00 TOL • •1000000000QE»1# 

XtU* .32331666552E 00 Xt2)- 18329745696E 01 


* INTERMEDIATE VALUES 16 GOOD STEPS TAKEN , ' 1 BAD STEPS TAKEN 

T • .30000000000E 01 FPT - .lOOOOOOOOOOE 01 DTP • .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 TOL • *10000000000E*1# 
X(ll« -.18660739124E 01 XC2)« • . 10210603358E 01 


INTERMEDIATE VALUES 14 GOOD STEPS TAKEN » 0 BAD STEPS TAKEN 

T ■ .40000000000E Ol FPT • .lOOOOOOOOOOE 01 DTP • .lOOOOOOOOOOE 01 DTQ T .lOOOOOOOOOOE 00 TBL • • lOOOOOOOOOOE* it 

XtU* ••17417683244E 01 XC2) * .624*66l*401E 00 


INTERMEDIATE VALUES 9 GOOD STEPS TAKEN $ 0 BAD STEPS TAKEN 

T - .50000000000E 01 FPT » .lOOOOOOOOOOE 01 DTP ■ .lOOOOOOOOOOE 01 DTG ■ .lOOOOOOOOOOE 00 TBL « .10000000000C*l# 

X cl) ■ ••837Q7745035E 00 XC2)* • 13070889378E 01 


INTERMEDIATE VALUES 15 GOOD STEPS TAKEN » 0 BAD STEPS TAKEN 

T • .60000000000E 01 FPT • .lOOOOOOOOOOE 01 DTP • .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 T0L • .100000000001*1# 

X(l3 * .1279Q420293E Ol X C2) ■ .24378144489E 01 


INTERMEDIATE VALUES 19 GOOD STEPS TAKEN # 0 BAD STEPS TAKEN 

t * .70000000000E 01 FPT • .lOOOOOOOOOOE 01 DTP • .lOOOOOOOOOOE 01 DTG • .lOOOOOOOOOOE 00 TBL * • 10000000000S*4# 
XCD* . 19201 524 166E 01 XC2) • ».435f 3I53332E 00 


C0 

CO 



CD 


TABLE A -4, (Continued) 


INTERMEDIATE VALUES 10 GOOD STEPS TAKEN . 0 BAD STEPS TAKEN 

T • .80000000000E Ot FPT « ,10000000000E 01 DTP • .iOOOOOOOOOOE 01 DTG • • 10000000000E 00 T0L • . lOOOOOOOOOOC»iO 

*03- • 12132324410E 01 X C2J • - , 98781392226E OQ 


INTERMEDIATE VALUES U GOOD STEPS TAKEN , 1 BAD STEPS TAKEN 

t • .90000000000E 01 PPT • .IOOOOOOOOOOE 01 DTP • .IOOOOOOOOOOE 01 DTG • .IOOOOOOOOOOE 00 TBL • • IOOOOOOOOOOE*10 

X Cl 3- -.A1291605205E 00 X C23 • • *25269034480E Ol 


INTERMEDIATE VALUES 17 GOOD STEPS TAKEN # 0 BAD STEPS TAKEN 

T - .IOOOOOOOOOOE 02 FPT • .IOOOOOOOOOOE 01 DTP • .IOOOOOOOOOOE 01 DTG • .IOOOOOOOOOOE 00 TBL *' • 1000000000QE«iO 

XClJ • -.20083407830E Ol XC2j« . 3290707Q692E-01 

INTERMEDIATE VALUES 13 GOOD STEPS TAKEN p 1 BAD STEPS TAKEN 

T • .llOOOOOOOOOE 02 FPT • .IOOOOOOOOOOE 01 DTP • .iOOOOOOOOOOE 01 DTG f .IOOOOOOOOOOE 00 TBL ■ • lOOOOOOOOOQC*lO 

Xll)» -.1S0A9739797E Ol Xt2}« ,78M 4M2429E 00 

INTERMEDIATE VALUES 13 G08D STEPS TAKEN # 0 BAD STEPS TAKEN 

T • .12000000000E 02 FPT • ♦ IOOOOOOOOOOE 01 DTP • . iOOOOOOOOOOE 01 DTQ • • IOOOOOOOOOOE 00 TBL ■ #1000000000qt*lS 

XCU ■ -.31376909536E 00 XC2]« .18440269613E Ol 


INTERMEDIATE VALUES 17 GOOD STEPS TAK£N , 0 BAO STEPS TAKEN 

T • .13000000000E 02 FPT • .IOOOOOOOOOOE 01 DTP • .IOOOOOOOOOOE 01 DTG • .IOOOOOOOOOOE 00 TBL • •lOOOOOOOQOOC.tO 

XtU ■ .18717708928E Ol XC2]« .99758656955E 00 


INTERMEDIATE VALUES 20 GBBD STEPS TAKEN , 0 BAD STEPS TAKEN 

T • .IOOOOOOOOOOE 02 FPT - .IOOOOOOOOOOE 01 DTP ■ .IOOOOOOOOOOE 01 DTG • .IOOOOOOOOOOE 00 TBL • . lOOOOOOOOOOf.tO 

K{l}m .1738641703GE Ol XC23- - .62715364635E 00 

INTERMEDIATE VALUES 8 GBBD STEPS TAKEN » 0 BAD STEPS TAKEN 

T • .15000000000E 02 FPT. .IOOOOOOOOOOE 01 DTP- • lOOOOOQOOOOE 01 DTG • • IOOOOOOOOOOE 00 TBL'* • 10000000000t*iO 



TABLE A -4. (Concluded) 


xtn- 

• 8304»37*2862E 

00 

XC23- 

- * 131336588A2E 

01 










intermediate VALUES 

U 

GOOD 

STEPS TAKEN , 0 

BAD STEPS TAKEN 

T ■ 

.Uoooooqoooe 

02 

FPT « 

•10000000000E 

01 OTP • .lOCOOOOOOOOE 

01 

DTG 

¥ .10000000000E .00 

, T9l ■ .1OOOOO0OOO0E.10 

xtn- 

-.12913899229E 

01 

X C2] ■ 

-.2423263800AE' 

01 











INTERMEDIATE values 

17 

GOOD STEPS TAKEN # 

1 BAD STEPS TAKEN 

T • 

•17000000C00E 

02 

FPT . 

.10000000000E 

01 DTP ■ • 1 OO n OOOOOOOE 

01 

DTG » .10000000000E 

00 TOl • tlOOOOOOOOOQE.lO 

xtn* 

-.19179362020E 

01 

X C2] • 

•A3962156AQ5E 

00 










INTERMEDIATE VALUES 

10 GOOD STEPS TAKEN , 

0 BAD STEPS TAKEN 

T • 

.18000000000E 

02 

FPT . 

.iooooooooooe: 

01 DTP ■ .lOCOOOOOOOOE 

01 DTG ■ .IOOOOOOOOOOE 

00 TOL • .1000000000QE*10 

xtn* 

-.12082152063E 

01 

X C 23 * 

.9316153i»7*2E 

00 









intermediate values 

10 

GOOD 

STEPS TAKEN , 0 BAD STEPS TAKEN 

T . 

.19000000000E 

02 

FPT • 

.iooooooooooe 

01 

DTP • 

.lOCOOOOOOOOE 

01 

DTG 

• .IOOOOOOOOOOE oo 

T8U • 

•lOOOOOOOOOOE.iO 

XClJ- 

•A257A876675E 

00 

XC2] • 

.253535A22A6E 

01 













FINAL 

VALUES 276G00D 

STEPS TAKEN 6BAD STEPS TAKEN 


T • 

.20000000000E 

02 

FPT « 

.iooooooooooe 

01 

dtp ■ 

.IOOOOOOOOOOE 

01 

DTG 

• .IOOOOOOOOOOE 00 

TBt • 

. iooooooooooe. io 

X(l]* 

•20081A9761AE 

01 

X C2] • 

• ♦ A2508886890E 1 

• 01 






... 




APPENDIX B 


A LISTING OF THE RATIONAL FUNCTION 
EXTRAPOLATION ALGORITHM 


A listing of both the old and new versions of DIFSYF is included here. 
These subroutines are versions of the rational function extrapolation technique 
as developed by Bulirsch and Stoer [5] . The listings are included here so that 
they may be compared with other versions that exist. 
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I 


Old Version of DIFSYF 


1 . 


SUBROUTINE DIFSYF (N, F, EFS, H,X, Y) 

2. 


DIMENSION Y( 17), DTT( 17) , D(6), YA(l7), YL(l7), YM(l7), 

3. 


DY(l7), D1(17), DT(l7,6), YG(8,17), YH(8,17), S(l7) 

4. 


INTEGER R, SR 

5. 


LOGICAL KONV, BO, BH, FIN 

6. 


ep=abs(eps) 

7. 


N1=N 

8. 


HH=H 

9. 


IF (EP.LT. 5. E-8) EP=5. E-8 

10. 


CALL F (X, Y, DZ) 

11. 


BH=. FALSE. 

12. 


FIN=. FALSE. 

13. 


DO 1 1=1, N1 

14. 


S(I)=0. 

15. 

1 

YA(I)=Y(I) 

16. 

2 

A=HH+X 

17. 


FC=1. 5 

18. 


BO=. FALSE. 

19. 


M=1 

20. 


R=2 

21. 


SR=3 

22. 


JJ=0 

23. 


DO 23 J=l,10 

24. 


IF (BO) GO TO 3 

25. 


D( l)=2. 25 

26. 


D(3)=9. 

27. 


D(5)=36. 

28. 


GO TO 4 

29. 

3 

D(l)=l. 7777777778 

30. 


D(3)=7. 1111111111 

31. 


D(5)=28. 4444444444 

32. 

4 

KONV=J. GT. 5 

33. 


IF (J. LE.7) GO TO 5 

34. 


L=6 

35. 


D(6)=64. 

36. 


FC=.6*FC 

37. 


GO TO 7 

38. 

5 

L=J-1 

39. 


IF (J-1) 7,7,6 

40. 

6 

d(l)=m*m 

41. 

7 

M=M+M 

42. 


g=hh/float(m) 
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43. 


B=G+G 

44. 


IF (BH.AND.J.LT.9) GO TO 14 

45. 


KK=(M-2)/2 

46. 


M=M-1 

47. 


DO 8 1=1, N1 

48. 


YL(I)=YA(I) 

49. 

8 

YM(l)=G*DZ(l)+YA(l) 

50. 


DO 13 K=1,M 

CALL F (X+FLOAT(K)*G,YM,DY) 

51. 


52. 


IF (DY(l).GT. 1.E38) GO TO 25 

53. 


DO 10 1=1, N1 

54. 


U=B*DY(l)+YL(l) 

55. 


yl(i)=ym(i) 

56. 


ym(i)=u 

57. 


U=ABS(U) 

58. 


IF (U-S(I)) 10,10,9 

59. 

9 

S(I)=U 

60. 

10 

CONTINUE 

61. 


IF (K.EQ.KK.AND.K.NE.2) GO TO 11 

62. 


GO TO 13 

63. 

11 

JJ=1+JJ 

64. 


DO 12 1=1, N1 

65. 


yh( jj,i)=ym(i) 

66. 

12 

yg(jj,i)=yl(i) 

67. 

13 

CONTINUE 

68. 


GO TO 16 

69. 

14 

DO 15 1=1, N1 

70. 


YM(I)=YH(J,I) 

71. 

15 

yl(i)=yg( J,l) 

72. 

16 

CALL F (A, YM,DY) 

73. 


IF (DY(l).GT. 1.E38) GO TO 25 

74. 


DO 22 1=1, N1 

75. 


V=DTT(I) 

76. 


DTT(I)=( YM(I) +YL(I) -K3*DY(l) )* . 5 

77. 


C=DTT(I) 

78. 


TA=C 

79. 


IF (L.LT.l) GO TO 20 

80. 


DO 19 K=l, L 

81. 


B1=D(K)*V 

82. 


B=B1-C 

83. 


u=v 
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84. 


IF (ABS(B)-l.E-lO) 18,18,17 

85. 

17 

B= ( C -V) /B 

86. 


U=C*B 

• 

00 


C=B1*B 

88. 

18 

v=dt(i,k) 

89. 


dt(i,k)=u 

90. 

19 

TA=U+TA 

91. 

20 

IF (ABS(Y(l)-TA).GT.EP*S(l)) GO TO 21 

92. 


GO TO 22 

93. 

21 

KONV=. FALSE. 

94. 

22 

Y(I)=TA 

95. 


IF ( KONV ) GO TO 24 

96. 


D(2)=4. 

97. 


D(4)=16. 

98. 


BO=.NOT.BO 

99. 


M=R 

100. 


R=SR 

101. 

23 

SR=M+M 

102. 


BH=.NOT. BH 

103. 


FIN=. TRUE. 

104. 


HH=HH*. 5 

105. 


GO TO 2 

106. 

24 

H=FC*HH 

107. 


X=A 

108. 


RETURN 

109. 

25 

H=0 

110. 


DO 26 1=1, N1 

111. 

26 

Y(I)=YA(I) 

112. 


END 


END OF COMPILATION: 


NO DIAGNOSTICS 



New Version of DIFSYF 


PrOTOkOLL 

subroutine 0IFSYS(N,F,EPS/H/X>Y) 
external f 

DIMENSION Y<17)*YA(21>*Yl<21>/YM<2n>0Y<21>/ 
1DZ<21),DT(2W)/D(7)j>S<21)*6PU> 

LOGICAL KO'IV*BU*GR*KL 

DATA EP/0.4E-l*0*16E-2*0.64E-A*D.256E-5/ 
JTI»0 

FY«1. 

ETAsABS < EPS ) 

IFIETA.LT.l.E-o) ETA-l.E-9 
DU IDO 1 = 1* 'I 
IOO YA ( I ) sy ( I ) 

CAUL F<X*Y*0Z) 

10 XN*X+H 

BU=. FALSE. 

DO llfi I = 1 .» I 

110 S(I)=C. 

, 1=1 
JR = 2 
JS = 3 

DH 260 J = l*l-> 

1F1.4UT.U0) GO TO 200 
»(2). l. 777777777 F8 
) (A) =7.11111111111 
L)(6) = 28 .444V*A«t4A 
Ui) T 1 ?0i 

200 D ( 2 ) = 2 . 2 5 
D ( A ) a 9 . 

U(6)=3t. 

201 1F1J.UF.7) r.,l TO 202 
L*7 

u ( 7 ) =64 . 

\ 00 TO 203 

202 L = J 
D(L)=M+M 

203 Knr.V = L.CT.3 
M = fWN 

0 = M/FLt.AT('l) 

B = G+G 

DU 2 1C 1 = 1 /N 
YL ( I ) =Y A ( I ) 

210 YM ( I ) = V A ( I 1 +G*OZ t I ) 

M =1-1-1 

DO 220 K = 1 * 

CALL F(X+FLnAT(«>*G,YM,DY) 

DO 220 I = 1*N 
U=YL ( I )+B*0Y( I ) 

YL< I )=VM( 1 1 
YM( I )«L 
U»ABS(L ) 

220 IF(U.CT.Sm) S(I)=U 
CALL F(XNjYH.jy) 

KL*L.LT.2 

0R=L.GT.5 

FS-O. 

DIJ 233 I = l*’l 
V«OT( I^D 

C«(Y'1< I ) +YL ( I ) +G*OY ( I ) >*0.5 
DT ( I * 1 ) *C 
TA = C 

IF(KL) GOTO 233 
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S02l*-r.J3 FKEUND 


DICKM2 OPF DFVUK-RZU 


U Li 23l K = 2 > L 
)*V 

)>=ai-c 
.*i = c — \/ 
i = V 

IFtB.EC.O. ) &L1TQ 230 
H = .1 / 3 
U = C*B 

t=ai«e 

*30 V»OT(I,K) 

"TU,!<)=U 
231 V tt s U + T A 

lH,JnT.K(' : 'V) TU 232 

IF { A3S (V ( 1 )-TA ) .GT.SU ) *ET A ) KUNV* * FALSE • 
232 IMi;r,i‘R.S{I).C0*0.) GO TO 233 

• V=AdS(k) /$< U 
IMFS.LT.f7) F.S-FV 
233 Y ( I ) = T A 

IMFS.I'Q.P. ) GOTO 250 

FA = FY 

*-L-l 

FY* ( Ep ( ► )/ FS )**< 1. /FLOAT (L+K) ) 

1 F ( L . E L . * ) G-lTil 240 
IF(FV.LT.0.7*FA) GO to 250 
24 . 1F(FY.GT.0.‘ 7 ) GOTO 2*0 
H=H#F Y 

jti=jti+i 

IF(JTI.GT. 3) GO TO 30 
GO T ) 10 

25u IF (KIM) GO TH 20 
U(3)=4. 

0(5) =16. 
y|J 3 . \it:i . B»j 
0 = JR 
JR = JS 

26o JS=H+M 
H = '1*0.5 
G-1 T’J 10 

.20 X = X ) 

H = II*FY 
GFTir*!, 

30 (1 = 0, 

I'D 3 00 I = 1 > i 
300- YU)=YA(I) 

RETURN 

LNU 
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